4

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 enter image description here

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.

enter image description here

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;
  • 1
    I don't think your projection is 3857, which is why there is no intersection. I have looked at the accompanying tif.xml file which suggests this is a file in Montana State Plane NAD83, which according to spatial ref sys is 2256, spatialreference.org/ref/epsg/2256. However, importing the file using -s 2256 also doesn't work. I have also opened the file in QGIS and looked at it, and under properties, the projection is stated as user defined, 100002, wth paramters +proj=lcc +lat_1=45 +lat_2=49 +lat_0=44.25 +lon_0=-109.5 +x_0=600000 +y_0=0 +datum=NAD83 +units=m +no_defs – John Powell Nov 30 '17 at 9:17
  • 1
    I have no experience worked with state plane data, but it seems clear that your issue is with the projection being used. The data is fine and there is no error in the SQL itself. – John Powell Nov 30 '17 at 9:18
  • 1
    Have you tried gdalwarp with the -s_srs epsg:2256 switch explicitly set. The problem is definitely with the projection and not the SQL. QGIS couldn't read the projection, so I suspect gdal can't either. – John Powell Nov 30 '17 at 16:51
  • 1
    So, no response to my question of Nov 30th about s_srs, which seemed quite reasonable. I have downloaded the file, checked the SQL, made comments, been ignored, and now the bounty (which frankly wasn't the motivation from the start) has vanished. I do sometimes wonder why I bother. – John Powell Dec 7 '17 at 17:42
  • 1
    Ok, sorry, no worries. I hope things improve for you. I have never seen a bounty vanish before -- I think you have broken the system :-) I am more than happy to try for you if you give me a link to the original file in 2256. – John Powell Dec 7 '17 at 21:56
3
+100

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
  • Does this not work in some way? It is always nice to get feedback after spending time to help someone. – John Powell Nov 8 '17 at 8:30
  • Hey John - This did not work. I am updating my question now to give you some more to work with. – jotamon Nov 28 '17 at 17:24
  • 1
    Hey, I think from what you say, something like: SELECT (ST_SummaryStats(ST_Clip(rast, geom, true))).* FROM tiledclassification64_2017-04-30, clip or whatever your clip geometry table is called, ie, do the clip inside the SummaryStats call. – John Powell Feb 1 '18 at 16:35
  • I am so close on this. Check out my update with the results from a clip derived st_summarystats – jotamon Feb 5 '18 at 15:44
  • 1
    I am glad this worked for you. Postgis raster is great, but it certainly takes some getting used to. I am looking forward to a beer in the Rockies one day :D, ideally after skiing. – John Powell Feb 7 '18 at 17:04

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.