I'm trying to calculate an average value in a buffer around a point.
Given a raster in EPSG 4326 and a similar point, this is fine, except that the values will be skewed as they are not projected in an equal-areas projection.
So - I've added to my spatial_ref_sys the SRID '97099' from saptialreference.org: http://spatialreference.org/ref/sr-org/7099/
This works fine for transforming points [ edit - added example ]:
-- get all the diff areas of diff pix values in one tile in which a point lies SELECT ST_Area((ST_DumpAsPolygons(r.rast)).geom), (ST_DumpAsPolygons(r.rast)).val, r.rid FROM rasters AS r INNER JOIN ( SELECT ST_Transform(ST_SetSRID(ST_MakePoint(-2, 51), 4326), 97099) as pt ) AS p ON ST_Intersects(r.rast, p.pt) AND r.filename = 'someraster.tif';
But using it with the raster fails with an error saying it cannot find SRID 97099 in spatial_ref_sys: [ edit - added error message ]
ERROR: rtpg_getSR: Cannot find SRID (97099) in spatial_ref_sys SQL state: XX000 Context: SQL function "st_transform" statement 1
SELECT ... WHERE SRID = 97099 confirms it is there... So I'm a bit lost. Is it querying from somewhere else, or maybe ignoring anything that isn't EPSG? I've followed the example here, for rasters:
Here's the full test SQL, the commented out lines show what I want to do, and the line above them will actually work (but is wrong, as it will still be 4326)
DROP TABLE IF EXISTS tmp_points; CREATE TEMP TABLE tmp_points ("gid" int, "point" geometry); INSERT INTO tmp_points ("gid", "point") VALUES (100231, ST_MakePoint(172.3, -43.76667)), (100233, ST_MakePoint(171.71667, -43.91667)), (100246, ST_MakePoint(170.92009, -44.09834)), (100290, ST_MakePoint(170.95009, -44.44835)), (100497, ST_MakePoint(173.28271, -42.77409)), (100498, ST_MakePoint(171.90983, -43.9236)), (1, ST_MakePoint(-2.136154, 49.220288)) ; -- get the average value arount a point for a given radius -- SELECT gid, filename, lon, lat, CAST(SUM(ST_Area((foo.gv).geom)*(foo.gv).val)/SUM(ST_Area((foo.gv).geom)) AS decimal(9,7)) AS avgimr, COUNT(foo.gv) as parts FROM ( SELECT ST_X(p.point)::NUMERIC(9, 5) AS lon, ST_Y(p.point)::NUMERIC(9, 5) AS lat, p.gid, sn.filename, ST_Intersection( sn.rast, -- THIS COMMENTED LINE BELOW FAILS, THE ABOVE WORKS (is wrong) -- ST_Transform(sn.rast, 97099, 'Bilinear'), ST_Buffer(ST_Transform(ST_SetSRID(p.point, 4326), 97099), 25000) ) AS gv FROM rasters sn, tmp_points p WHERE ST_Intersects( ST_Buffer(ST_Transform(ST_SetSRID(p.point, 4326), 97099), 25000), sn.rast -- THIS COMMENTED LINE BELOW FAILS, THE ABOVE WORKS (is wrong) -- ST_Transform(sn.rast, 97099, 'Bilinear') ) AND sn.filename = 'someraster.tif' ) AS foo WHERE (foo.gv).val >= 0 AND (foo.gv).val <= 10 GROUP BY filename, gid, lon, lat ORDER BY filename, gid, lon, lat;