1

I have several GeoTIFFs (.tif) that I am conducting raster calculations on using gdal_calc.py, all three are of different sizes. These rasters can be best illustrated with the following image:

Raster Domains

All three rasters domains overlap across the northern extent, however the blue raster is of full size and covers the entire area of interest.

I'd like to perform a simple raster calculation (say gdal_calc.py -A redgrid.tif -B greengrid.tif -C bluegrid.tif --calc="(A+B+C)/3"), and generate an output raster in which I will get resultant values across the extent of the largest raster available, or the blue domain. Currently, this calculation will only yield an output raster with the smallest extent, or the red domain.

I've considered filling the other two grids with np.nan to match the extent of the blue domain, though I believe that this division would return nan in areas where the three grids do not overlap.

8
  • gdal-translate all layers to the same extent, resolution, CRS and origin first. That will put nodata values in empty cells. Make sure nodata values are defined. Then try the raster calculation. All from the command line.
    – wingnut
    Commented Apr 28, 2021 at 17:55
  • @wingnut Thanks for looking into this! These 3 rasters actually went through this process prior to the final calculation, all have an identical extent, XY resolution, CRS, & origin to each other, with a NoDataValue set to -9999. Is that where my problem exists at this time? Commented Apr 28, 2021 at 18:06
  • The result is normal. It only works on valid data. Change your equation for cases where only 1 or 2 layers exist.
    – wingnut
    Commented Apr 28, 2021 at 18:10
  • Ok good to know! I assume using a numpy conditional wrapped in gdal_calc.py {using an example from other code I've written: divide(A,B,out=zeros_like(A), where=B!=0)} wouldn't work in this case either? Commented Apr 28, 2021 at 18:14
  • If the equation is (A+B+C)/3") why not to expand all images into same size and fill extra cells with 0 but without NoDataValue?
    – user30184
    Commented Apr 28, 2021 at 18:33

1 Answer 1

2

Zeroing will disturb the mean.

Maybe a more complex solution, case-by-case is called for. If nodata is -9999

Cases

There are 7 cases:

  • 1 case where all 3 exist
  • 3 cases where 1 doesn't but 2 do exist
  • 3 cases where 1 does but 2 don't exist

Test if all 3 exist. * acts like a logical AND if arguments are 1 or 0:

(A>-9999) * (B>-9999) * (C>-9999)

produces 1 or 0 for each pixel. Multiply that by (A+B+C):

(A>-9999) * (B>-9999) * (C>-9999) * (A+B+C)

to get the answer where all 3 pixels are defined. If only 2 pixels are defined, the third must be OFF.

(A>-9999) * (B>-9999) * (C=9999) * (A+B)

There are 2 more cases like this. When only one dataset is valid at a point, the same applies

(A>-9999) * (B=-9999) * (C=-9999) * A

There are 2 more cases like this.

Combining Cases

Just as * represents AND, use + to represent OR. There seem to be 7 cases (add as required). Final equation is thus

(A>-9999) * (B>-9999) * (C>-9999) * (A+B+C) +
(A>-9999) * (B>-9999) * (C=9999) * (A+B) +
(A>-9999) * (B=-9999) * (C>9999) * (A+C) +
(A=-9999) * (B>-9999) * (C>9999) * (B+C) +
(A>-9999) * (B=-9999) * (C=-9999) * A +
(A=-9999) * (B>-9999) * (C=-9999) * B +
(A=-9999) * (B=-9999) * (C>-9999) * C

Line breaks are just to keep it neater. Maybe use a text editor to write it and then paste it into gdal_calc_py :-)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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