5

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?

  • Could you clarify "the original geometry intact"... – Pierre Racine Sep 12 '18 at 12:46
  • Is your original raster table tiled? – Pierre Racine Sep 12 '18 at 12:49
  • @PierreRacine Original geometry intact means that the rasterized output of the ST_MapAlgebra averages into a band must contain the same point geometries relative to the source SRS of the original raster. Essentially, the center pixels of the output band must match the geographic location of the pixels in the source band. I was a bit worried that in the process of using the ST_MapAlgebra function and unionizing the band I would end up with a no-geometry raster. – Evan McCoy Sep 12 '18 at 14:30
  • @PierreRacine The raster is tiled at 100x100 on import from the GeoTiff – Evan McCoy Sep 12 '18 at 14:31
  • @PierreRacine. What difference does tiling make to the mechanics of MapAlgebra. I assumed it was more of a performance issue. – John Powell Sep 12 '18 at 16:31
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.

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.