2

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: example

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.

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.

| improve this answer | |
1

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.

| improve this answer | |
1

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.

| improve this answer | |
  • 1
    Option 1 won't work - 'Dump your raster as polygons'. There would be 33 billion polygons. I've tried this with points and get memory errors (even with 64GB of RAM). Polygons would be worse due to more complex geometry. Thanks for suggestion 2, I'll look into it! – Theo F Jan 30 at 15:24
  • Ok. then use the second one. The difficulty of solution 2 is creating the callback function. You have to pass all the parameters of a raster to your callback function. Then, inside the callback function you will create that raster again just to have the position of each pixel. then, you will use that pixel (can be in geom format) to do your calculations. The only example I found online is the link that I sent you. You have to pass all the parameters that the code in line 1179 passed to the function. – milad Jan 30 at 15:48

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.