4

I have several raster tables in PostGIS (four in fact but in this case let us assume there are just two) - same cell sizes, same alignment but with different types of value (e.g 3d_pot, 3d_area). Occasionally there will be -9999 values where data is missing but this might not be the case in all three datasets for the same cell (e.g. there could be 3d_pot but not 3d_area).

I want to calculate the aggregate statistics for all the cells that are within a polygon (e.g. a building footprint) but only for those cells where 3d_pot is greater than a value (say 100).

For one value I need to calculate a combination of the cell values at the same location (3d_area * 3d_pot).

I'm struggling to find a fast way to effectively join the rasters based upon their spatial location, having intersected with the polygon layer and ignore all cells if one of them is -9999.

What is the simplest SQL way of doing this on the fly (without generating new rasters or rasterizing the building polygons?)

Here is what I currently have but this feels wrong and excessive to do multiple Intersections. Can you join one raster to another just using the raster geometry?

WITH solarpot AS (
SELECT fid, 
    (ST_Intersection(T.wkb_geometry, ST_Clip(ST_MapAlgebra(AREA.rast,POT.rast,'([rast2])',NULL,'INTERSECTION'), 1, T.wkb_geometry, True))).val AS potval,
    (ST_Intersection(T.wkb_geometry, ST_Clip(ST_MapAlgebra(AREA.rast,POT.rast,'([rast1])',NULL,'INTERSECTION'), 1, T.wkb_geometry, True))).val AS areaval,
    (ST_Intersection(T.wkb_geometry, ST_Clip(ST_MapAlgebra(AREA.rast,POT.rast,'([rast1]*[rast2])',NULL,'INTERSECTION'), 1, T.wkb_geometry, True))).val AS sumpotval
   FROM solar_rasters.solarp_3d_area AREA
   JOIN os_gla.topographicarea T ON ST_Intersects(T.wkb_geometry, AREA.rast)
   JOIN solar_rasters.solarp_3d_pot POT ON ST_Intersects(T.wkb_geometry, POT.rast)
   WHERE T.fid = 'osgb5000005156542071'
)
SELECT fid, COUNT(potval), SUM(potval), SUM(areaval) AS sumareaval, SUM(sumpotval) AS sumpotval, AVG(potval), stddev(potval), MIN(potval), MAX(potval)
FROM solarpot
WHERE potval IS NOT NULL AND potval >= 100 AND areaval IS NOT NULL
GROUP BY fid
  • Are your rasters tiled? – Pierre Racine Oct 4 '19 at 13:09
  • Yes the are: raster2pgsql -s 27700 -t 100x100 -I -C -M – SpatialSDE Oct 6 '19 at 14:27
0

So I have currently done this by just making a JOIN on the raster geometry and this seems to work reasonably well (given the size of the rasters). If I am missing something basic and there is a much faster way to do this it would be good to know:

-- this works and seems to get us correct values
WITH solarpot AS (
   SELECT fid, 
    (ST_Intersection(T.wkb_geometry, ST_Clip(POT.rast, 1, T.wkb_geometry, True))).* as geom

   FROM solar_rasters.solarp_3d_pot POT
   JOIN os_gla.topographicarea T ON ST_Intersects(T.wkb_geometry, POT.rast)
   WHERE T.fid = 'osgb1000005305136'
),
area AS (
   SELECT fid, 
    (ST_Intersection(T.wkb_geometry, ST_Clip(AREA.rast, 1, T.wkb_geometry, True))).* as geom

   FROM solar_rasters.solarp_3d_area AREA
   JOIN os_gla.topographicarea T ON ST_Intersects(T.wkb_geometry, AREA.rast)
   WHERE T.fid = 'osgb1000005305136'
),
tilt AS (
   SELECT fid, 
    (ST_Intersection(T.wkb_geometry, ST_Clip(TILT.rast, 1, T.wkb_geometry, True))).* as geom

   FROM solar_rasters.solarp_3d_tilt TILT
   JOIN os_gla.topographicarea T ON ST_Intersects(T.wkb_geometry, TILT.rast)
   WHERE T.fid = 'osgb1000005305136'
),
azi AS (
   SELECT fid, 
    (ST_Intersection(T.wkb_geometry, ST_Clip(AZI.rast, 1, T.wkb_geometry, True))).* as geom

   FROM solar_rasters.solarp_3d_azi AZI
   JOIN os_gla.topographicarea T ON ST_Intersects(T.wkb_geometry, AZI.rast)
   WHERE T.fid = 'osgb1000005305136'
)
SELECT SP.fid, SP.val AS potval, A.val AS area3d, TI.val AS tilt, AZ.val AS azi
FROM solarpot SP 
JOIN area A ON SP.geom = A.geom
JOIN tilt TI ON SP.geom = TI.geom
JOIN azi AZ ON SP.geom = AZ.geom
WHERE SP.val > 500

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.