Tuesday, 25 July 2017

geotiff tiff - Use OSGeo4W shell command from python 2.7 script


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

Answer



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.



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...