1

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

  • If you are using GDAL >= 2.1, use the gdal.Translate() function which avoids having to call a subprocess. – user2856 Mar 13 '18 at 21:47
  • I tried: gdal.Translate(tif_file,jp2_file) but I get ValueError: Received a NULL pointer. – user88484 Mar 14 '18 at 7:31
1

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.

  • 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? – user88484 Mar 16 '18 at 9:22
  • 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. – jdbcode Mar 16 '18 at 15:34
1

@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.

0

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.

  • how does it work? should I write : 'gdal_translate -a_srs ESPG:32639 -of Gtiff input.jp2 output.tif' where the callstring is located? – user88484 Mar 13 '18 at 21:24
  • @user88484 Sorry, thought that was obvious. But my first time through, nothing was obvious :P See edit. – Jon Mar 13 '18 at 21:28
0
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 ESPG:4326 -of Gtiff input.jp2 output.tif')

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

  • 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. – user88484 Mar 13 '18 at 21:22
  • @user88484 Try the update and let me know if it works. – Marcelo Villa Mar 13 '18 at 23:06
  • 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 – user88484 Mar 14 '18 at 7:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.