3

I am looking to run the equivalent of Focal Statistics in gdal_calc. Is it possible?

It is a simple focal I require. I need to run a 3 by 3 kernel over the data (mean) and if possible only on NoData (but if on all that is okay as well).

I can do this in other software but specifically asking if a GDAL_calc solution exists.

This is the exact thing I am trying to do on a singular raster

Con(IsNull("C:\temp\SST_Final_C_NoData.tiff"), FocalStatistics("C:\temp\SST_Final_C_NoData.tiff", NbrRectangle(5,5, "CELL"), "MEAN"), "C:\temp\SST_Final_C_NoData.tiff")

GDAL 1.11.0 and Python 3.4.1

  • I have not tried myself but it should be possible to use kernel filters inside GDAL virtual rasters. Read gdal.org/gdal_vrttut.html and try. – user30184 Jun 10 '14 at 18:14
  • I think a python way exists through calc but cannot track it down. – If you do not know- just GIS Jun 10 '14 at 20:31
  • Gdal_calc is mostly undocumented except in code, see trac.osgeo.org/gdal/ticket/5388. I hope you can find a way to do what you want. You can also ask from gdal-dev mailing list which is also for this kind of questions. – user30184 Jun 10 '14 at 20:56
4

I don't believe you can do it directly through gdal_calc unfortunately (though I'm happy to be proved wrong). However you can use the gdal_calc.py script as a basis to read the data into Python and do the operation through scipy.ndimage which has a pile of inbuilt focal filter operations. Unfortunately it doesn't have a mean filter, but if you're happy with median you could use:

result = scipy.ndimage(your_raster_as_numpy_array, size=3)

Alternately if you do need the mean you could use a convolution filter and divide the result by the number of cells covered:

kernel = numpy.ones((3,3))
result = scipy.ndimage.convolve(your_raster_as_numpy_array, weights=kernel) / kernel.size

You should then be able to write that back out to a file in the same way the gdal_calc script does.

A couple of notes though:

  • Edge effects are still going to be an issue. Happily the filter methods in scipy.ndimage have a number of ways to deal with these
  • Large rasters may be problematic - you might need to read them in by chunks and process them that way. You will have to deal with edge overlaps though, so that you avoid even more edge effects
  • No data gaps are likely to be a pain. You will have to think about how you're going to deal with holes in your raster, either by interpolating to fill the gaps, or by setting the nodata values to numpy.NaN, which you will then have to interpret when writing back out through gdal

Good luck, and hope it helps!

1

I think I have found it or close enough.

gdal_fillnodata SST_Final_C_NoData.tiff -md 5 -si 2 SST_Final_C_NoData_Foc.tiff

Fills No Data with a 5 pixel window (-md 5) and does 2 passes (-si 2)

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.