Get ST_SummaryStats from Raster Table using geoJSON polygon

Question

I have a raster dataset imported into POSTGIS with SRS 3857. This table is named as prcp0000_normal3857.

The dataset represents statewide precipitation in Montana and can be downloaded below. ftp://mco.cfc.umt.edu/prcp/Daymet/yearly_normals/geotiff/PRCP0000_normal.zip

After downloading the data was reprojected to 3857, tiled and renamed prcp0000_normal3857 using the GDAL command below.

gdalwarp -co "TILED=yes" -of GTiff -t_srs epsg:3857 -r bilinear prcp0000_normal.tif prcp0000_normal3857.tif

The data was loaded into the POSTGIS database test1 using the command below

raster2pgsql -I -C -s 3857 /dataFolder/prcp0000_normal3857.tif -F -t 100x100 | psql -d test1

I also have a geojson polygon in SRS 4326 that I would like to use to get some zonal statistics on the raster

The geojson looks like this

'{"type":"Polygon","coordinates":
             [[[-108.797607421875,47.73193447949174],
             [-109.31396484375,47.234489635299184],
             [-107.9296875,46.800059446787316],
             [-106.89697265625,47.12247581664114],
             [-107.02880859375,47.76148371616669],
             [-108.797607421875,47.73193447949174]]]

}'

Based on suggestions from John's answer I am trying the POSTGIS command below

WITH clip(geom) AS 
  (SELECT 
     ST_Transform(
       ST_SetSRID(
         ST_GeomFromGeoJSON(
       '{"type":"Polygon","coordinates":
           [[[-108.797607421875,47.73193447949174],
                 [-109.31396484375,47.234489635299184],
                 [-107.9296875,46.800059446787316],
                 [-106.89697265625,47.12247581664114],
                 [-107.02880859375,47.76148371616669],
                 [-108.797607421875,47.73193447949174]]]
             }'),
        4326), 
     3857)
   )
 SELECT ST_Clip(rast, geom, true) 
   FROM prcp0000_normal3857, clip;

POSTGIS just returns the errors below with no other info

I feel like I am close but not indicating the band and perhaps messing up the ST_Transform command

Update: I can now successfully run the ST_Summary stats on the table and make the clip. I just need to know how to bring two together.

The ST_SummaryStats quickly returns the mean on the entire table by using the comment below. I just need to know how to pass the clipped raster piece instead of the entire "prcp0000_normal3857" table

SELECT (ST_SummaryStats(rast)).mean FROM "prcp0000_normal3857";

Some more updates. I've fixed projection issues etc and can successfully run ST_SummaryStats using the intersect between the raster and the geojson. That looks like this

  SELECT (ST_SummaryStats(rast)).*
  FROM "prcp0000_normal3857"
  WHERE st_intersects(rast, (ST_Transform(
    ST_SetSRID(
      ST_GeomFromGeoJSON(
        '{"type":"Polygon","coordinates":
            [[[-116.8402862548828,37.59464778787345],[-116.75582885742188,37.59464778787345],[-116.75582885742188,37.67050562029702],[-116.8402862548828,37.67050562029702],[-116.8402862548828,37.59464778787345]]]
          }'),
     4326),
  3857)));

This returns summary stats below but not exactly clipped since I used st_intersects. Results are below. I am thinking there are multiple rows because my raster is tiled.

Bringing it all back together. When I try and clip the raster and run st_summarystats, I still get an error that the raster and polygon do not intersect. At this point, the data are correctly projected and I feel I am missing something in the syntax of this command.

WITH clip(geom) AS
  (SELECT
     (ST_Transform(
       ST_SetSRID(
         ST_GeomFromGeoJSON(
           '{"type":"Polygon","coordinates":
               [[[-116.8402862548828,37.59464778787345],[-116.75582885742188,37.59464778787345],[-116.75582885742188,37.67050562029702],[-116.8402862548828,37.67050562029702],[-116.8402862548828,37.59464778787345]]]
             }'),
        4326),
     3857))
   )
SELECT (ST_SummaryStats(ST_Clip(rast, geom, true))).* FROM 
"prcp0000_normal3857", clip;

Answer 1

One thing I often do is check the formats available for a particular function, which you can do in psql using:

\df+ ST_Clip

which returns among other things, an argument data type list, one of which is:

rast raster, geom geometry, crop boolean

which does not involve specifying any bands. The 2nd example in ST_SummaryStats docs

ST_SummaryStats(ST_Clip(rast,2,geom))

appears to be out of date -- the only version I can see that matches that would be (rast, band [], geom, nodata, crop), ie, which an array of bands.

So, returning to your example, and moving the geometry out into a CTE to make it easier to see what is going on:

WITH clip(geom) AS 
  (SELECT 
     ST_Transform(
       ST_SetSRID(
         ST_GeomFromGeoJSON(
           '{"type":"Polygon","coordinates":
               [[[-107.07275390625001,40.613952441166596],
                 [-107.66601562500001,40.195659093364654],
                 [-107.22656250000001,39.53793974517628],
                 [-106.10595703125001,40.1452892956766],
                 [-105.66650390625001,40.697299008636755],
                 [-107.07275390625001,40.613952441166596]]]
             }'),
        4326), 
     3857)
   )

 SELECT ST_Clip(rast, geom, true) 
   FROM testdataedited, clip;

works for me, albeit on a local dataset I have.

You have an extra pair of parenthesis before ST_Transform, which is not an error as such, but detracts from readability a little. The error:

NOTICE: The input raster and input geometry do not intersect. Returning empty raster

might suggest other issues, but also might go away once you have the correct form of ST_Clip.

EDIT: Following from the comments and having downloaded the source Tiff and opened it in QGIS -- which showed unknown projection -- I believe the issue lies in the not setting the source CRS in the gdal translate command.

gdalwarp -co "TILED=yes" -of GTiff -t_srs epsg:3857 -r bilinear prcp0000_normal.tif prcp0000_normal3857.tif

should be:

gdalwarp -co "TILED=yes" -of GTiff -s_srs:epgs:2256 -t_srs epsg:3857 -r bilinear prcp0000_normal.tif prcp0000_normal3857.tif

EDIT 2 following exchange in comments and by mail.

The final working query is:

WITH clip(geom) AS
  (SELECT
    (ST_Transform(
      ST_SetSRID(
        ST_GeomFromGeoJSON(
          '{"type":"Polygon","coordinates":
             [[[-116.8402862548828,37.59464778787345],
               [-116.75582885742188,37.59464778787345],
               [-116.75582885742188,37.67050562029702],
               [-116.8402862548828,37.67050562029702],
               [-116.8402862548828,37.59464778787345]]]
           }'),
      4326),
   3857))
 )
SELECT (ST_SummaryStats(ST_Clip(rast, geom, true))).* 
  FROM "elevtestfloattiled", clip 
 WHERE ST_Intersects(rast, geom);

The important addition here is the inclusion of the final WHERE clause. Because the input raster is tiled, it is necessary to use ST_Intersects to only select those tiles that intersect the clip -- otherwise you get lots of errors about non-overlapping geometries. Personally, I think this is one of the confusing things about raster coverages: it seems like one table, with one set of constraints, but, via rid, is actually made up of many separate rows, which constitute one logical coverage.

The results from this look like:

count |       sum        |       mean       |      stddev      |       min        |       max        
  230 | 391537.942526309 | 1702.33888054917 | 43.8034115158831 | 1622.91694047492 | 1806.35017687067
  506 | 979477.957614989 | 1935.72718896243 | 109.492359004463 | 1725.76388786231 | 2237.12976290935
  130 | 216469.197384467 | 1665.14767218821 |  30.085111261015 | 1616.02460735199 | 1722.30092577995
  286 | 534027.315403459 | 1867.22837553657 | 90.8440344578266 | 1718.70971908873 | 2078.68385757086

Finally, if you want all the clipped statistics in one line, you can wrap the ST_Clip in ST_Union, ie, change the last line to:

SELECT (ST_SummaryStats(ST_Union(ST_Clip(rast, geom, true)))).* 
  FROM "elevtestfloattiled", clip
 WHERE ST_Intersects(rast, geom);

which now returns:

 count |       sum        |       mean       |      stddev      |       min        |       max        
  1152 | 2121512.41292922 | 1841.59063622328 | 138.708602370447 | 1616.02460735199 | 2237.12976290935