Intersection of a MultipolygonZ and a floating point raster

Question

I am trying to isolate the pixels in a floating point raster (loaded in PostGIS) that intersect my MultipolygonZ table (also loaded in PostGIS) and whose values are less than or equal to the z values of the MultipolygonZ table - essentially a 3d intersection.

I am running PostGIS 2.1.5. SELECT postgis_full_version() returns:

"POSTGIS="2.1.5 r13152" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER"

Here's what I have tried so far:

1) ST_Intersects(rast, geom) - this simply calculates a 2d intersection between the polygons and the raster - the z values of the MultipolygonZ are ignored.

SELECT (ST_DumpAsPolygons(rast)).*
FROM pgz_tab, dem
WHERE ST_Intersects(rast, geom)

2) ST_3dIntersects(rast, geom) - invalid arguments - apparently this only works with points and linestrings - http://postgis.net/docs/manual-2.1/ST_3DIntersects.html

SELECT (ST_DumpAsPolygons(rast)).*
FROM pgz_tab, dem
WHERE ST_3dIntersects(rast, geom)

3) Converting the MultipolygonZ to a floating point raster and then subtracting the two rasters. I can't seem to get the arguments right with the ST_AsRaster function - http://postgis.net/docs/manual-2.1/RT_ST_AsRaster.html even after reading A guide to the rasterization of vector coverages in PostGIS.

Answer 1

You thought this through well, and it might be that raster subtraction is your best bet. But, just to leave it registered, as of PostGIS 2.2 ST_3dIntersects() does work with polygons and TINs. Creating said polygon3d or TIN from a raster in PostGIS is another problem entirely, but once it's done, the function will work just fine.

Regarding your problem specifically: there are many ways to work with ST_AsRaster(). Different parameters you can call, etc. In your case, the easiest might be the second-to-last one on the list, in which you only need to set 4 parameters:

ST_AsRaster(geom, width, height, pixel_type)

where geom is your geometry object, width and height the raster's dimensions, and pizel_type the PostGIS' Raster Data Type. Beware, however, that this will create only an empty raster, in which every pixel is set to value 1. You have to then set the values.

For that, you can use the ST_SetValues() function. It takes in a raster object, the raster band (first band is 1), and a geomval[] type. GeomVal is an aggregate function that holds a set of geometries, and a set of values of each geometry - you can think of it as a point cloud, though you can construct it with linestrings and polygons as well. You construct it with Potsgres' array_agg (you can check an implementation of it here). Something like:

SELECT ST_SetValues(ST_AsRaster('POLYGON((0 0, 5 0, 5 5, 0 5, 0 0))'::geometry, 5, 5, '8BUI'::text), 1, geomval)
FROM (SELECT array_agg((p.geom, p.val)::geomval) as geomval
      FROM public.poly3d as p
) values

Quite a workload, to be frank, but if you can't upgrade your system to PostGIS 2.2 or above, you may have to make do with this.