GDAL VRT Pixel Function no output with signed or float data type

Question

I am trying to average multiple partially overlapping rasters using a VRT and a Python based Pixel Function (https://gdal.org/drivers/raster/vrt.html#using-derived-bands-with-pixel-functions-in-python). The rasters have cells with no data.

The function on this page was my starting point: https://gis.stackexchange.com/a/324680/169329. The problem I am running into is that I am unable to change the dataType to one that allows me to have fractions and negative values without having useless output.

I have tried multiplying by 1000 and rounding the values before building the VRT and averaging the values using a Pixel Function, but this works only for positive values with dataType UInt32.

My expectation is that it should be possible to use other dataTypes in Pixel Functions besides those that only contain positive and rounded numbers. How can I do this?

I have used the following code after first multiplying the rasters by 1000 to try and get an average before dividing by 1000 again:

<VRTRasterBand dataType="Int32" band="1" subClass="VRTDerivedRasterBand">
<PixelFunctionType>average</PixelFunctionType>
<PixelFunctionLanguage>Python</PixelFunctionLanguage>
<PixelFunctionCode><![CDATA[
    import numpy as np

    def average(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,raster_ysize, buf_radius, gt, **kwargs):
        div = np.zeros(in_ar[0].shape, dtype = 'Int32')
        for i in range(len(in_ar)):
            div += (in_ar[i] != 0)

        div[div == 0] = 1
        y = np.sum(in_ar, axis = 0, dtype = 'Int32')
        y = y / div
        out_ar[:] = np.round(np.clip(y,-100000,100000))
]]>
</PixelFunctionCode>

Answer 1

While I see this is an old post, I recently had the exact same issue and solved it like this:

<VRTRasterBand dataType="Float32" band="1" subClass="VRTDerivedRasterBand">
  <PixelFunctionType>average</PixelFunctionType>
  <PixelFunctionLanguage>Python</PixelFunctionLanguage>
  <PixelFunctionCode><![CDATA[
    import numpy as np

    def average(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,raster_ysize, buf_radius, gt, **kwargs):
        x = np.ma.masked_equal(in_ar, -3.4e+38)
        np.nanmean(x, axis = 0,out = out_ar, dtype = 'float64')
        mask = np.all(x.mask,axis = 0)
        out_ar[mask]=-3.4e+38
    ]]>
  </PixelFunctionCode> 

The code first masks any NoData. The no data value of -3.4e+38 is the default for a Float64 GDAL raster and this was my clue as to why it wasn't working - that and the overflow warning message. Then, using np.nanmean allows me to do the calculation in a single step without first counting the number of overlaps. I finally replace the masked values with my no data value.

However, I have to explicitly set dtype to 'float64' and when I used GDAL translate to 'fix' the result, I also had to explicitly define the input type also as Float64. Any variation from these settings resulted in the overflow warning and nonsense in the output.