Proper use of python subprocess for GDAL

Question

I need to process some Sentinel-2 files which are in JP2 format using python. I first want to clip the JP2 files and save the result as GeoTIFF with the same CRS. But some of the modules I use (e.g., rasterio for clip) don't like jp2 format so I try to convert the jp2 to GeoTIFF. I am able to do so using QGIS translate tool which is essentially GDAL translate. However, when I try to use the translate tool in Python:

import geojson, gdal, subprocess
tif_file = r"C:\Users\Administrator\L2A_B04_10m.tif"
jp2_file = r"C:\Users\Administrator\L2A_T29RNP_20170422T110651_B04_10m.jp2"
args = ['gdal_translate', '-of', 'Gtiff', jp2_file , tif_file]
subprocess.Popen(args)

I get a file without CRS. I tried to embed the CRS in the args but I can't figure out a way to do so, I don't get any error but it seems that the code is ignoring my tries, like:

args = ['gdal_translate', '-a_srs','ESPG:4326', '-of', 'Gtiff', jp2_file , tif_file]
subprocess.Popen(args)

results: nothing happens, no error, and no file is created

also:

args = ['gdal_translate', '-a_srs ESPG:4326', '-of', 'Gtiff', jp2_file , tif_file]
subprocess.Popen(args)

also:

args = ['gdal_translate', '-a_srs', '+proj=longlat +ellps=WGS84', '-of', 'Gtiff', jp2_file , tif_file]
subprocess.Popen(args)

and:

args = ['gdal_translate', '-a_srs +proj=longlat +ellps=WGS84', '-of', 'Gtiff', jp2_file , tif_file]
subprocess.Popen(args)

What am I missing? or, is there another easier way to clip JP2 format using Python?

I'm using Python 2.7 on anaconda on Win10, GDAL 2.2.3

• user2856: If you are using GDAL >= 2.1, use the gdal.Translate() function which avoids having to call a subprocess.
• user88484: I tried: gdal.Translate(tif_file,jp2_file) but I get ValueError: Received a NULL pointer.

Answer 1

import os
import subprocess
os.chdir("working directory")
# Option 1 List of arguments
args = args = ['gdal_translate', '-a_srs','ESPG:4326', '-of', 'Gtiff', jp2_file, tif_file]
result = subprocess.call(args)
# Option 2 String
result = subprocess.call('gdal_translate -a_srs ESPG:4326 -of Gtiff input.jp2 output.tif')

You can either pass a string or a list of arguments. However try calling subprocess.call() instead of subprocess.Popen(). In the string, you'll have to change the file names respectively.

Update:

Try calling gdal_warp as follows:

subprocess.call('gdal_warp -t_srs EPSG:4326 -of Gtiff input.jp2 output.tif')

Additionally, there are other creation options -co you might be interested in. Check the following link: http://www.gdal.org/gdalwarp.html

• user88484: both options resulted in nothing.. no error, no tif file. When I delete the coordinate part 'ESPG:4326', I do get a file but without coordinate system which causes problems later.
• Marcelo Villa: @user88484 Try the update and let me know if it works.
• user88484: does not work.. it has the same results as the gdal_translate. It seems that the problem is not the specific gdal tools, but rather the use of the srs in the syntax. I can't find the right way to use it

Answer 2

@user88484 - I would suggest you get the commands working at the command prompt (just use one sample file) before you try to automate/batch it. I do this to disentangle issues with the commands from automation errors.

I do something similar in a two-step process. First, I register the data to the coordinate space I want by providing the bounding box. In this case, I am reading a NetCDF to a Tif.

 gdal_translate -b 5 -of GTiff -a_ullr <ULX ULY LRX LRY> -a_nodata 0 -co COMPRESS=LZW NETCDF:in.nc:PM25 out.tif

Then, I use gdalwarp to reproject it to another projection (it's actually to a geographic coordinate system with a different spheroid if you want to get technical). Here, I define both the input and the output.

gdalwarp -s_srs +proj=longlat +a=6371200 +b=6371200 +no_defs -t_srs EPSG:4326 -srcnodata "value 0" -overwrite in.tif out.tif

Note that in other cases, I have resorted to using the -co PROFILE=BASELINE option with gdal_translate. My understanding is that this totally wipes out all the Tif tags (projection and otherwise) and just gives you a Tif with the values and you have to build it back up. This is a scorched-earth approach, but drastic times... See the profile section about 2/3 down the GDAL GTiff page.

As for subprocess call vs. Popen, see here. Note that GDAL produces a lot of warnings in the STDERR. Look at the returncode from subprocess more than the STDERR to determine success/failure.

Answer 3

When I get a file written without CRS using a subprocess call to a GDAL util, the problem is usually solved by explicitly stating where the proj4 directory can be found using a config setting, as in --config GDAL_DATA <path to */Anaconda2/Library/share/gdal>, mine is in: C:/Users/user/Anaconda2/Library/share/gdal. Try something like:

args = ['gdal_translate', '-of', 'Gtiff', '--config GDAL_DATA C:/mock/path/Anaconda2/Library/share/gdal', jp2_file , tif_file]
subprocess.Popen(args)

You can also set the GDAL_DATA path as an environmental variable.

• user88484: I tried it and nothing happens, no error, no new file. My folder is C:\ProgramData\Anaconda2\Library\share\gdal I also tried to search for some folder with proj4 and found C:\ProgramData\Anaconda2\pkgs\proj4-4.9.3-vc9_5 but again, it resulted in nothing. Any other suggestion?
• Justin Braaten: Bummer. It would be helpful to know more about what is returned from your call to gdal_translate. subprocess won't return everything by default, so try what @Jon suggested or run the command in command prompt. These should provide better warning and error messaging.

Answer 4

If you want to see your errors (and outputs), use

args = ['gdal_translate', '-a_srs','ESPG:4326', '-of', 'Gtiff', jp2_file , tif_file]
proc = subprocess.Popen(args, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
stdout,stderr=proc.communicate()

stderr contains error messages, stdout contains output messages.

Also note that "-a_srs" overrides whatever projection gdal reads by default from your file. It does not reproject your output raster. Use gdalwarp for that.

• user88484: how does it work? should I write : 'gdal_translate -a_srs ESPG:32639 -of Gtiff input.jp2 output.tif' where the callstring is located?
• Jon: @user88484 Sorry, thought that was obvious. But my first time through, nothing was obvious :P See edit.