Populate an empty raster with values from intersecting vector layers

Question

I'm looking to populate a large, empty Unsigned 8-bit raster with integer values 1 to 10. The value depends on the spatial relationship between a given raster cell and a variety of vector layers (points/lines/polygons) stored in PostGIS. I also have the vectors in shapefile and File GeoDatabase formats. Various fields and attributes within the vector layers need to be taken into account too.

Example: IF raster cell is 50m away from polygon A AND polygon A name = 'Alpha', THEN assign value of 3 to that raster cell. And repeat for all cells etc....

Bearing in mind the empty raster is 5x5m and covers the entire area of Great Britain (so billions of cells), and that the vector layers (30-40 or so) can contains 1000s of features each, what is the best way to populate the raster with values?

Raster details: -ot Byte -of GTiff -tr 5 5 -a_nodata 255

Example of the process I'm looking to repeat:

Note how this isn't simply a case of rasterizing vector layers, it's about the spatial relationship between multiple layers as well. See how the red dots, although small, actually affect a greater area proportionally than say, the blue polygon.

Answer 1

I would recommand to process the spatial relationships within the vector format, then rasterize the result. It seems that some buffers around your vectors combined with priority rules would work quite well. All the operations that I know between raster and polygon in fact imply a format conversion at one stage. So either you create the buffers first, then rasterize and apply some rules. Or you rasterize, then create the buffers with a dilation filter and apply your rules.

Answer 2

It should be possible to do with gdal_rasterize https://gdal.org/programs/gdal_rasterize.html#gdal-rasterize.

You could for example create a set of buffered geometries for your source data. Burn little values with the largest geometries first, and then increase the pixel values by burning the smaller geometries into the raster by using the -add option

-add

Instead of burning a new value, this adds the new value to the existing raster. Suitable for heatmaps for instance.

Answer 3

I can think of two solutions:

1. Dump your raster as polygons using this function ST_PixelAsPolygons(). Then, do your calculations for the centroid of each polygon (first calculate the centroids of the polygons). Now, you should have a table of points that each of them has a value that you calculated. Finally, build an array of your points (a geomval object) and use this function ST_SetValues() to assign the values of the points to your empty raster.

2. The second solution is to use this function ST_MapAlgebra(). This function runs for each pixel of your raster and allows to pass a function to it and it will run that function for each pixel and assign the output of the function to that pixel. If you want to see an example of that function you can see this code https://github.com/pedrogit/postgisaddons/blob/master/postgis_addons.sql line 1179 a callback function is created.

The second solution seems to be faster (I think) but creating the callback function may be demanding. Let me know if you want more help for creating the callback function.