Use OSGeo4W shell command from python 2.7 script

Question

I would like to convert jp2 files (sentinel-2) to geotiff format using python 2.7 on win 10. I saw some solutions using python and GDAL but they did not work for me. The only thing that did work was to open OSGeo4W shell, go to the desired directory using cd directory (from the shell) and then run (in the shell) for %i in (*jp2) do gdal_translate -of GTiff -co TFW=YES %i %~ni_tiff.tif. The problem is that I need to perform it on multiple folders and also there are some processes I need to perform on each file before and after the conversion, the conversion from jp2 to geotiff is just a small part of the process. So I would like to know if it is possible to call OSGeo4W shell command from within a python script. My desired results would be a python script that at some part of the script runs a command on the OSGeo shell, each time for a different folder (I can store the folders in a list if necessary) and continues to process the products of the OSGeo4W command shell.

I read solutions regarding command line (cmd) using os.system or subprocees but none of them referred to executing commands in a different platform than win command line.

The script according to my logic would be something like that:

import ...
#src is a list of directories that holds directories with jp2 files
src = ["c:\jp2_files_1", "c:\jp2_files_2","c:\jp2_files_3","c:\jp2_files_4"]
for i in (src):
 os.chdir(i)
 SomethingI'mMIssing('for %i in (*jp2) do gdal_translate -of GTiff -co TFW=YES 
 %i %~ni_tiff.tif')
continue with python processing

• Dmitry Baryshnikov: It seems to me you need to start you script from OSGeo4W shell, not start shell from script.

Answer 1

Adding to what @bennos said, I needed to run some commands in Windows PowerShell and ran into a similar problem. You can easily switch to OSGeo4W-Shell by inserting its program call to your actual call, like this:

import subprocess as sub
# adjust path if necessary
my_call = [r'C:\Program Files\QGIS 2.18\OSGeo4W.bat',
 'gdalbuildvrt',
 '-input_file_list',
 vrt_textfile,
 vrt_file]
# call it
p = sub.Popen(my_call, stdout=sub.PIPE, stderr=sub.PIPE)
stdout, stderr = p.communicate()
if p.returncode != 0:
 print stdout
 print stderr

---

To be a bit more specific regarding your actual problem (using gdal_translate in a loop on many files), you could do it like this:

import os
import subprocess as sub
# get your files
path = r'c:\MY_PATH'
files = [os.path.join(path, f) for f in os.listdir(path) if f.endswith('.jp2')]
# now run gdal_translate on each file
for f in files:
 basename = os.path.splitext(os.path.basename(f))[0]
 newname = os.path.join(path, basename + '_translated.tif')
 # set up call
 my_call = [r'C:\Program Files\QGIS 2.18\OSGeo4W.bat',
 'gdal_translate', '-of', 'GTiff',
 '-co', 'TFW=YES',
 f, newname]
 # call it
 print 'Using gdal_translate on {0} to convert it to {1}'.format(f, newname)
 p = sub.Popen(my_call, stdout=sub.PIPE, stderr=sub.PIPE)
 stdout, stderr = p.communicate()
 if p.returncode != 0:
 print stdout
 print stderr

Instead of using this way of getting your filenames you could also just create a fixed list of files as you did initially. subprocess (or sub in my case) will open the OSGEO4W Shell in the background and perform the given my_call.

I assign this to the variable p to be able to retrieve its return code (0 if it executed successfully) to see if something went wrong and print not only the error message (stderr), but also the "normal" messages produced by gdal_translate (stdout).

You can test the script by leaving out the actual call (everything below # call it) and use print my_call instead just to see what the actual would look like.

• user88484: Hi, I'm not an expert on the issue so I would really appreciate if you can elaborate more, more step-by-step. I didn't understand if I can run the script @s6hebern suggested in the python platform (jupyter notebook), also, it seems to me that you did not use the my_call variable in your script, and is the cmd is a variable or something else? thanks
• s6hebern: @user88484 - See edits. I'm not sure about the possibility to run this in Jupyter, since never worked with it, but just give it a try.
• user88484: Amazing.. thank you @s6hebern. At first I got error message regarding if p.exitcode != 0:. The error was: 'Popen' object has no attribute 'exitcode'. But since the actual code works I deleted these "communicating" lines
• s6hebern: My bad, it actually is p.returncode. Edited my answer. Glad it works :)

Answer 2

I use it exactly like @DmitryBaryshnikov commented: calling GDAL as a subprocess in the python script and then executing the pyhon-file in the OSGeo4W shell. It is like your Win command line.

Some example code:

my_call = [
 'gdalbuildvrt',
 '-input_file_list',
 vrt_textfile,
 vrt_file, 
]
subprocess.call(my_call)