GDAL calculation without stop python.exe error?

Question

I want to create a simple raster calculation using python and GDAL and specific gdal_calc.py in osgeo4w shell.

BUT any time that before finish my script python stopping I take error : python.exe stop working.

I have read to many post from users with this problem but without specific solution.

in my PATHEXT I missed '.PY' I try to add but without success.

here the first code :

gdal_calc =  'python C:/OSGeo4W64/bin/gdal_calc.py -A path/to/a.tif -B path/to/b.tif  --outfile=path/to/out.tif --calc="((4*A)+(10*B))+4.8"'

try to add pathext :

call("set PATHEXT = % PATHEXT %;.PY;.PYW")
gdal_calc =  'python C:/OSGeo4W64/bin/gdal_calc.py -A path/to/a.tif -B   path/to/b.tif  --outfile=path/to/out.tif --calc="((4*A)+(10*B))+4.8"'

in this case I get windows error 2

after I try this :

os.environ['PATHEXT'] = '.PY'
gdal_calc =  'python C:/OSGeo4W64/bin/gdal_calc.py -A path/to/a.tif -B   path/to/b.tif  --outfile=path/to/out.tif --calc="((4*A)+(10*B))+4.8"'

but in this case gdal_calc.py don't execute .

Any idea how to can work this script without error ?

Answer 1

Here is a template I use when calling gdal tools from python:

import subprocess, sys
callstr = [sys.executable,
           'C:/OSGeo4W64/bin/gdal_calc.py', 
           '-A', 'path/to/a.tif',
           '-B', 'path/to/b.tif',
           '--outfile=path/to/out.tif'
           '--calc=((4*A)+(10*B))+4.8'\
           ]
proc = subprocess.Popen(callstr, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
stdout,stderr=proc.communicate()

Some things to note: When using subprocess, you shouldn't put quotes around your calc equation [I am not sure this is true; read it on another answer I can't find now]. The output variable stdout will contain any strings returned by the process, and stderr will return any errors. Because gdal_calc.py isn't an executable, you need to specify its path and include the sys.executable command, apparently.

Finally, if your rasters are relatively small, depending on your processing needs you might find it more useful to load them as numpy arrays and do your raster math with numpy. If this route sounds interesting, look at rasterio which makes reading/writing geotiffs pretty easy.