Generating a neighborhood pixel raster from existing raster using PostGIS

Question

I am starting with a PostGIS raster table containing anywhere from 6-12 bands. In a nutshell, the goal is to take each of these bands, perform a 3x3 neighborhood average on each pixel in the band, and output a new 12-band raster with the neighborhood-averaged bands and the original geometry intact. I could likely do this with pure GDAL, but the rasters at this stage in the process have already been imported into the database and the original GeoTIFF files deleted, so using GDAL would likely require an infrastructure change.

My current concept uses ST_MapAlgebra to take the bands from the source raster in one-by-one and perform the 3x3 average to return the single averaged band. I then would use ST_Union on a setof these averaged bands to create a single raster. I then need to somehow generate a new table from the output of the ST_Union. I am trying to do this in as few steps as possible, so ideally the process would take 1-2 queries.

My ST_MapAlgebra callback function is this:

CREATE OR REPLACE FUNCTION
                    callback_fn(pixel float[], pos integer[], variadic userargs text[])
                        returns float
                        language plpgsql
                        immutable
                    as $$
                        declare
                            avg float;
                            x integer;
                            y integer;
                            p integer[];
                        begin
                            x = pos[1][1];
                            y = pos[1][2];
                            p = pixel[1];
                            avg := (p[x-1][y-1] + p[x][y-1] + p[x+1][y-1]
                                    + p[x-1][y] + p[x][y] + p[x+1][y]
                                    + p[x-1][y+1] + p[x][y+1] + p[x+1][y+1]) / 9.0;

                            return avg;
                        end
                    $$;

My current attempt at an execution query is definitely not syntactically correct, but it gets the point across of what I'm trying to do:

WITH rast AS (
                        SELECT rast FROM "raster2"
                    )
                    CREATE TABLE "raster3" as (select ST_Union(
                        ARRAY[st_mapalgebra(rast, 1, 'callback_fn(float[], integer[], text[])'::regprocedure),
                        st_mapalgebra(rast, 2, 'callback_fn(float[], integer[], text[])'::regprocedure)]
                    ));

It is likely I am way over-complicating the process, but I am fairly new to PostGIS. One of my concerns is that an approach like this will not preserve the underlying geometry of the raster, which is a must. How can I change my queries to be functional in the way I need. Aside from that, am I even approaching this in the right way, or is there a much cleaner way?

Answer 1

The answer I came up with didn't end up using processing tools within PostGIS, but instead used GDAL. It turned out to be a fairly efficient method, even though it had to convert to an intermediate GeoTIFF file to function.

raster = gdal.Open("""PG:""" + config.connstring() + """ schema='public' table='table_name' column='rast' mode='2'""")
outFile = gdal.GetDriverByName("GTiff").Create("neighborhoodRaster.tif", raster.RasterXSize, raster.RasterYSize, raster.RasterCount*2, gdal.GDT_Float32)
outFile.SetGeoTransform(raster.GetGeoTransform())
outFile.SetProjection(raster.GetProjection())

for bandId in range(1, raster.RasterCount):
    band = raster.GetRasterBand(bandId).ReadAsArray()
    outFile.GetRasterBand(bandId).WriteArray(band)
    outFile.GetRasterBand(bandId + raster.RasterCount).WriteArray(neighborhoodStatistics(band))

def neighborhoodStatistics(data):
        data = np.pad(np.array(data, dtype=np.float32), pad_width=1, mode='constant', constant_values=np.nan)
        shape = (data.shape[0] - 2, data.shape[1] - 2) + (3, 3)
        strides = data.strides + data.strides
        windows = np.lib.stride_tricks.as_strided(data, shape=shape, strides=strides)
        return np.nanmean(windows, axis=(2, 3))

This solution reads a raster from the database and outputs a GeoTIFF file with the same georeference, but double the bands. The second set of bands is the 3x3 neighborhood mean of the raw set.

The trick to this was using Numpy and its built in stride_tricks module. The original raster is padded with a ring of NaN values to allow the edge pixels to have a full 3x3 scope. Stride_tricks is used to restride the band data array to a four-dimensional array of 3x3 neighborhood segments, one for each original pixel. A Numpy nanmean is then applied over the two axes representing the 3x3 neighborhood zones to reduce the band array back to two dimensions.

Using stride_tricks over a generic iterative window function made this process orders of magnitude faster! I was able to process a size band raster with ~7000x7000 pixels in a matter of seconds on my 8700k machine.