6

I try to assign each Polygon in my data set the maximum value that can be found within their areas. I have a raster and a collection of Polygons with each have a unique id (gid).

I uploaded the raster with:

raster2pgsql -d -I -C -M -F -t auto -s 32637 <Mypath> <Rastername> | psql -d <DBname> -U <Username>

When I run following code I get some results which appear to be reasonable for the segments, but I also get a lot of NULL values. When I run zonal statistics on the same geometry in QGIS (Raster Analysis > Zonal statistics) I get no NULL values while running GRASS v.rast.stats gives me exact the same output like in PostGIS, additional with the warning:

WARNING: Not all vector categories converted to raster. Converted 2128 of 3361.

Which correspondents also with the numbers of missing values in Postgis. I'm breaking my head whats causing this problem. I guess its related to the geometry created with ST_Buffer, wich somehow is not compatible with St_SummaryStats(). I also thought the problem is related to the tiling but even without tiling it occurs. ST_IsValid() returns for every element of my geometry TRUE, also is every element according to GeometryType() a Polygon.

  WITH blocks AS (
    SELECT row_number() OVER() AS gid, 
           ST_Buffer(ST_Line_Substring(the_geom, 8*n/length,
             CASE
               WHEN 8*(n+1) < length THEN 8*(n+1)/length
               ELSE 1
             END), 4, 'endcap=flat join=mitre') AS the_geom
    FROM (SELECT ST_LineMerge(geom) AS the_geom,
              ST_Length(geom) As length
              FROM <LineDataset>
         ) AS t
    CROSS JOIN generate_series(0,CEIL(length/8)::int) AS n
    WHERE n*8/length < 1
  )
  SELECT gid,
         (ST_SummaryStats(St_Union(ST_Clip(rast, 1, the_geom, true)))).max AS a
  FROM <Rastername>, blocks
  WHERE ST_Intersects(rast, the_geom)
  GROUP BY gid;

The following picture shows the attribute table of the output of QGIS zonal statistics (left) and the code from above (right).

enter image description here enter image description here

Some example data to recreate the error:

  1. Some random Landsat data
  2. Upload raster with: raster2pgsql -d -I -C -M -F -t 100x100 -s 3763 <somepath>\\Landsat8_L1TP_RGBN.tif test1 | psql -d <DBname> -U <Username>
  3. Create sample line e.g.:
CREATE TABLE line AS 
SELECT  
ST_SetSRID(ST_MakeLine(ST_MakePoint(-28718,172580),    ST_MakePoint(-14012,182541)), 3763)::GEOMETRY(LINESTRING, 3763) AS geom
  1. Run code from above after setting <Rastername> to test1 and <LineDataset> to line.
4
  • 1
    Could you provide a subset of your data or something similar? It may help for debugging the query.
    – wfgeo
    Commented May 30, 2020 at 17:50
  • Please have a look, I updated my question.
    – DrSnuggles
    Commented Jun 1, 2020 at 16:42
  • Do some of the polygons overlap?
    – Dataform
    Commented Jun 7, 2020 at 22:02
  • @Dataform I don't have time to help diagnose this, but the extents of the clipped rasters do overlap. Here's a screengrab.. The larger squares correspond to the envelopes of the raster segments with one each associated to one of the smaller squares.
    – Rob Skelly
    Commented Jun 7, 2020 at 22:12

1 Answer 1

7
+250

PostGIS rasters are more involved to work with than e.g. in QGIS, as QGIS does much of the work under the hood.

PostGIS stores raster pixels internally as grid references to the upper left pixel coordinate, and uses either the raster meta data (e.g. pixel dimension) in conjunction to work in 'raster space', or derives geometries from those values per pixel to work in 'vector space'.

Now, considering your example data, there are several issues evolving out of the fact that your polygons are a lot smaller than the pixel dimension of the raster:

  • ST_Intersects works on raster extents, and in the case of tile, on their individual extents; it is not finding intersections between geometries and pixels!

  • ST_Clip extracts 'pixels' by checking the input geometry against their centroids; since those buffer polygons are smaller than the pixel dimensions, only those pixels that happen to have their centroid falling inside a buffer polygon are extracted!

Thus, ST_SummaryStats is working only on those raster parts that were created from the few pixels that have been extracted by chance of their centroids falling within a buffer polygon.


When working with rasters in 'vector space', you should better fully think in terms of vector geometries, and vectorize the raster data accordingly; but keep in mind that vectorized rasters are very, very large, and the index on the raster column is not working, so your input geometries should come from a table with a spatial index in place! In every case, the process is going to take some time.

There are now several ways to get to your stats; As an example I'm going to use ST_PixelAsPolygons to dump pixels as a SETOF geomval, which holds the pixel envelope as polygon and the respective band value to aggregate over based on intersection with the buffer polygons; note that I assume your buffered polygon are in a table bp, with a spatial index in place:

WITH
  rclip AS (
    SELECT  ST_Clip(r.rast, 1, ST_Buffer(p.geom, ST_ScaleX(r.rast)), TRUE) AS rast
    FROM    test1 AS r
    JOIN    (
        SELECT  ST_Union(bp.geom) AS geom
        FROM    bp
    ) AS p  
      ON    ST_Intersects(r.rast, p.geom, 1)
  )

SELECT bp.gid,
       MAX(cells.val)
FROM   rclip AS rc,
       LATERAL ST_PixelAsPolygons(rc.rast, 1) AS cells
JOIN   bp
  ON   bp.geom && cells.geom
GROUP BY
       1
;

Here,

  • I first clip the raster to a buffer around the input geometries, with a radius corresponding to the X dimension of a pixel; this massively limits the amount of pixels to create polygons for (each tile in your example Landsat raster results in 10000 polygons). Unfortunately, to get a single raster in return, one has to either union the clipped rasters, or, as above, the input polygons into a single one. This would be faster if used with a buffer around the initial lines!
  • I then compute the maximum of all intersecting pixel polygons per input buffer polygon, grouped by their gid

This takes around 30 seconds on my machine, using your sample. I haven't done any performance hunting, and didn't compare other methods.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.