Calculate the derivative of a raster band

[IvanSanchez]

Feature description

Earlier this week (during FOSS4GE 2024 2025), I wanted to take a DEM and calculate a U-V slope map from it.

I know I can do gdal raster slope then gdal raster aspect then gdal raster stack then gdal raster calc to convert the slope-aspect polar coordinate into a u-slope,v-slope cartesian coordinate (or x-y slope). But that is doing trigonometry back and forth and feels like extra computing.

So I'd like to propose gdal raster derivative.

It should take a single band as input, and output two bands: the x-derivative and the y-derivative. The x-derivative band is, for each pixel, the difference between the pixels to the left and to the right; the y-derivative is the difference between the pixels up and below.

In other words: the x-derivative band is like applying this image kernel to the input band:

[ 0 0 0]
[-1 0 1]
[ 0 0 0]

...and the y-derivative band is like applying this kernel:

[0 -1 0]
[0  0 0]
[0  1 0]

Perhaps this tool could output a single band by having a parameter to discern between the X and Y directions. Otherwise I believe this would be a nice simple tool with its own use cases.

I've thought about specifying the angle for the derivative, or specifying custom kernels, but that sounds like overkill.

Additional context

This is just an idea, and I have no immediate need for it; please consider this feature request as low priority.

And greetings from the FOSS4G-EU codesprint!

[IvanSanchez]

Perhaps related: #11828

[msmitherdc]

@IvanSanchez can this be done via a pipeline and then stream to gdal raster calc and then output to GDALG?

[IvanSanchez]

can this be done via a pipeline and then stream to gdal raster calc and then output to GDALG?

@msmitherdc As far as I'm aware (and as far as I understand the documentation), the numPy expression in gdal raster calc cannot access neighbouring pixels.

I mean, I could clip one (or two) rows and columns, then use that and the original raster in gdal raster calc, but that would duplicate the temporary data.

[jratike80]

Can't you do that with kernel filtered VRT?
https://gdal.org/en/stable/drivers/raster/vrt.html#kernelfilteredsource

Use the same SourceBand, and create two output bands by applying different kernels. No duplication, no temporary data.

EDIT

I experimented a bit and the VRT might look about like this

<VRTDataset rasterXSize="3000" rasterYSize="3000">
  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["EUREF-FIN / TM35FIN(E,N)",GEOGCS["EUREF-FIN",DATUM["EUREF-FIN",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1391"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","10690"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",27],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3067"]]</SRS>
  <GeoTransform>  2.5400000000000000e+05,  2.0000000000000000e+00,  0.0000000000000000e+00,  6.9540000000000000e+06,  0.0000000000000000e+00, -2.0000000000000000e+00</GeoTransform>
  <VRTRasterBand dataType="Float32" band="1" blockXSize="256" blockYSize="256">
    <NoDataValue>-9999</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <KernelFilteredSource>
      <SourceFilename relativeToVRT="1">n3424h.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="3000" RasterYSize="3000" DataType="Float32" BlockXSize="256" BlockYSize="256" />
      <Kernel normalized="0">
        <Size>3</Size>
        <Coefs>0 0 0 -1 0 1 0 0 0</Coefs>
      </Kernel>
    </KernelFilteredSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="2" blockXSize="256" blockYSize="256">
    <NoDataValue>-9999</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <KernelFilteredSource>
      <SourceFilename relativeToVRT="1">n3424h.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="3000" RasterYSize="3000" DataType="Float32" BlockXSize="256" BlockYSize="256" />
      <Kernel normalized="0">
        <Size>3</Size>
        <Coefs>0 -1 0 0 0 0 0 1 0</Coefs>
      </Kernel>
    </KernelFilteredSource>
  </VRTRasterBand>
</VRTDataset>

The two derived bands look like this in QGIS, though slightly different. That is a desired result, I think.

I believe that the two bands can be combined into single band also with VRT by creating a VRTDerivedRasterBand https://gdal.org/en/stable/drivers/raster/vrt.html#vrt-derived-1. Specifying custom kernels is up to you.

Maybe you can do a prototype of gdal raster derivative tool yourself by modifying the gdalfilter.py sample script https://github.com/OSGeo/gdal/blob/master/swig/python/gdal-utils/osgeo_utils/samples/gdalfilter.py.

BTW, QGIS edits the VRT file when the layer is removed from the project, which is a bit irritating when testing. Adding the image statistics makes sense, but first I tried to use normalized coefficients -0.5 and 0.5, and for some reason QGIS always changes the "Kernel normalized" value into 0, without changing the coefficients.

[jratike80]

If there are MapServer folks in the FOSS4GE 2025 code sprint, they might be able to use such VRT as input for a layer of type "uvraster" https://mapserver.org/input/vector/vector_field.html and render slope arrows on-the-fly from DEM.

[rouault]

dding the image statistics makes sense, but first I tried to use normalized coefficients -0.5 and 0.5, and for some reason QGIS always changes the "Kernel normalized" value into 0, without changing the coefficients.

I cannot reproduce this. This would be a GDAL behavior anyway. What is the GDAL version of your QGIS ?