Build Mosaic using gdal based on Modal value

Question

I have an immense (~50,000) set of tiles that I want to mosaic using the mode for every pixel. The tiles have a size of 299x299 and are overlapping each other by 2/3 respectively.

For the sake of overview I only use 2 tiles now, which you can find here: https://drive.google.com/open?id=1A14hjemFtXhjIPib8amifTLmieNKm-G4

I tried mosaicking them using the Pixel-Function in the VRT as proposed in: Averaging overlapping rasters with gdal

So first I built a VRT with

gdalbuildvrt mosaic.vrt *.png

And in the Pixel-Function Code I tried it in two ways:

First:

        <VRTDataset rasterXSize="399" rasterYSize="299">
          <GeoTransform>  8.1042600000000000e+05,  3.0000000000000000e+00,  0.0000000000000000e+00,  2.4950340000000000e+06,  0.0000000000000000e+00, -3.0000000000000000e+00</GeoTransform>
          <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
            <PixelFunctionType>mode</PixelFunctionType>
            <PixelFunctionLanguage>Python</PixelFunctionLanguage>
            <PixelFunctionCode><![CDATA[
    import numpy as np
    from scipy import stats

    def mode(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,raster_ysize, buf_radius, gt, **kwargs):
        out_ar = stats.mode(in_ar, axis=0, nan_policy = 'omit')[0]
    ]]>
        </PixelFunctionCode> 
            <ColorInterp>Gray</ColorInterp>
            <SimpleSource>
              <SourceFilename relativeToVRT="1">prediction_92530.PNG</SourceFilename>
              <SourceBand>1</SourceBand>
              <SourceProperties RasterXSize="299" RasterYSize="299" DataType="Byte" BlockXSize="299" BlockYSize="1" />
              <SrcRect xOff="0" yOff="0" xSize="299" ySize="299" />
              <DstRect xOff="0" yOff="0" xSize="299" ySize="299" />
            </SimpleSource>
            <SimpleSource>
              <SourceFilename relativeToVRT="1">prediction_92531.PNG</SourceFilename>
              <SourceBand>1</SourceBand>
              <SourceProperties RasterXSize="299" RasterYSize="299" DataType="Byte" BlockXSize="299" BlockYSize="1" />
              <SrcRect xOff="0" yOff="0" xSize="299" ySize="299" />
              <DstRect xOff="100" yOff="0" xSize="299" ySize="299" />
            </SimpleSource>
          </VRTRasterBand>
        </VRTDataset>

Second:

    <VRTDataset rasterXSize="399" rasterYSize="299">
      <GeoTransform>  8.1042600000000000e+05,  3.0000000000000000e+00,  0.0000000000000000e+00,  2.4950340000000000e+06,  0.0000000000000000e+00, -3.0000000000000000e+00</GeoTransform>
      <VRTRasterBand dataType="Byte" 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):
        scores = np.unique(np.ravel(in_ar))
        oldmostfreq = np.zeros(in_ar[0].shape)
        oldcounts = np.zeros(in_ar[0].shape)
        for score in scores:
            template = (in_ar == score)
            counts = np.expand_dims(np.sum(template, 0),0)
            mostfrequent = np.where(counts > oldcounts, score, oldmostfreq)
            oldcounts = np.maximum(counts, oldcounts)
            oldmostfreq = mostfrequent
        out_ar = mostfrequent
    ]]>
        </PixelFunctionCode>    
        <ColorInterp>Gray</ColorInterp>
        <SimpleSource>
          <SourceFilename relativeToVRT="1">prediction_92530.PNG</SourceFilename>
          <SourceBand>1</SourceBand>
          <SourceProperties RasterXSize="299" RasterYSize="299" DataType="Byte" BlockXSize="299" BlockYSize="1" />
          <SrcRect xOff="0" yOff="0" xSize="299" ySize="299" />
          <DstRect xOff="0" yOff="0" xSize="299" ySize="299" />
        </SimpleSource>
        <SimpleSource>
          <SourceFilename relativeToVRT="1">prediction_92531.PNG</SourceFilename>
          <SourceBand>1</SourceBand>
          <SourceProperties RasterXSize="299" RasterYSize="299" DataType="Byte" BlockXSize="299" BlockYSize="1" />
          <SrcRect xOff="0" yOff="0" xSize="299" ySize="299" />
          <DstRect xOff="100" yOff="0" xSize="299" ySize="299" />
        </SimpleSource>
      </VRTRasterBand>
    </VRTDataset>

Then I took gdal_translate to calculate the mode as defined in the Pixel-Function Code:

gdal_translate --config GDAL_VRT_ENABLE_PYTHON YES mosaic.vrt raster.tif

However, the output is not correct. It is correct in size but it only returns me a raster full of zeros. Is there something I am overlooking or is there another way to mosaick the tiles based on their modal values? I tried both functions in Spyder and there I got the correct results, so I think that I am missing something at the end with the out_ar.

Answer 1

Here is the best solution I found up to now, however if someone has an idea on how to do it differentely or what I could do to improve the code I would appreciate it.

The Pixel-Function Code for calculating the modal value in the end is:

<VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
        <PixelFunctionType>average</PixelFunctionType>
        <PixelFunctionLanguage>Python</PixelFunctionLanguage>
        <PixelFunctionCode><![CDATA[
import numpy as np
from scipy import stats
def average(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,raster_ysize, buf_radius, gt, **kwargs):
    mode_arr = stats.mode(in_ar, axis=0, nan_policy = 'omit')[0][0]
    np.clip(mode_arr, 0,255, out = out_ar)
]]>
        </PixelFunctionCode>

However, one has to pay attention, as the numpy array reads the nodata values as real values. You would have to change them to NA first, and then to calculate the modal value.