1

I am trying to correct some water table .tifs (WT) to the ground level elevation (DEM) so that water levels are not incorrectly presented as above ground. I have looked at many similar questions and answers but it seems to me that the conditional capability of gdal_calc is limited to subbing in a static floating point number. I need it to use values from a DEM in the event that the water table is above ground. Pythonically, this would be

if (DEM - WT) > 0:
   output = WT
else:
   output = DEM

However, I can only seem to get gdal_calc to sub in a static floating point number, not a reference to a more dynamic .tif file where the value would depend on the grid location.

Is this possible with gdal? I have tried the following command line calls:

python "C:\Program Files\GDAL\gdal_calc.py" -A DEM.tif -B WT.tif --calc="B*((A-B)>0)" --outfile WT_belowGround.tif

This produces a .tif with values only where WT is below ground with No Data gaps where it is above ground. This is close to what I want but I want data coverage in all locations for producing contours.

python "C:\Program Files\GDAL\gdal_calc.py" -A DEM.tif -B WT.tif --calc="B*((A-B)>0)" --calc="A*((A-B)<0)" --outfile WT_aboveGround.tif

This produces a .tif with values only where WT is above ground with No Data gaps where it is below ground. Again, this is close to what I want but I want data coverage in all locations and it seems the second --calc call overwrites the first with No Data values.

Trying to combine both of these I have tried several things without the desired result. The most obvious thing to me is to assign the NoDataValue to the desired DEM.tif as is shown below.

python "C:\Program Files\GDAL\gdal_calc.py" -A DEM.tif -B WT.tif --calc="B*((A-B)>0)" NoDataValue=DEM.tif --outfile diff.tif

While this runs and produces a .tif file, it gives an error and again, it produces a .tif with values only where WT is below ground with No Data gaps where it is above ground. The error is below

RuntimeWarning: invalid value encountered in greater

Is this something that gdal can do or do I need to find another way to approach this? One thought I have is to produce the two .tif files of WT above ground and below ground and combine them but I am not sure how to combine them either.

3

In addition to expressions you can use numpy functions in gdal_calc. Sounds like you might be looking for np.where().

Something like:

--calc="where((A-B)>0), A, B)"
|improve this answer|||||
  • That is awesome! I was not aware of this, now I think I will be investigating numpy functions to see what other gems it contains. – Michael Stimson Feb 20 at 2:13

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.