# Generating a neighborhood pixel raster from existing raster using PostGIS

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"... Sep 12, 2018 at 12:46
• Is your original raster table tiled? Sep 12, 2018 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. Sep 12, 2018 at 14:30
• @PierreRacine The raster is tiled at 100x100 on import from the GeoTiff Sep 12, 2018 at 14:31
• @PierreRacine. What difference does tiling make to the mechanics of MapAlgebra. I assumed it was more of a performance issue. Sep 12, 2018 at 16:31

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):
outFile.GetRasterBand(bandId).WriteArray(band)
outFile.GetRasterBand(bandId + raster.RasterCount).WriteArray(neighborhoodStatistics(band))

def neighborhoodStatistics(data):