Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

A 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: http://www.postgis.org/documentation/manual-svn/RT_ST_Transform.html

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;

Many thanks

share|improve this question

1 Answer 1

up vote 1 down vote accepted

Looking at the Postgis 2.0.0 source, it seems the raster version of ST_Transform is only checking projections where the auth_name is 'EPSG'. See line ~512 in raster/rt_pg/rt_pg.c:

snprintf(sql, len, "SELECT CASE WHEN upper(auth_name) = 'EPSG' AND 
length(auth_srid::text) > 0 THEN upper(auth_name) || ':' || auth_srid ELSE '' END,
proj4text, srtext FROM spatial_ref_sys WHERE srid = %d LIMIT 1", srid);

See the source here: http://trac.osgeo.org/postgis/changeset/9787/trunk/raster/rt_pg/rt_pg.c

I'm not sure why it's limiting to EPSG, but it might be a bug worth filing on the project.

share|improve this answer
    
That is the correct reason, but in changing the auth_name and auth_srid to EPSG and 97099 respectively, the query now just blows up and disconnects from the DB without an error... :-( –  Steve Pike Jun 12 '12 at 14:06
    
For reference: trac.osgeo.org/postgis/ticket/1862 –  Steve Pike Jun 13 '12 at 10:06
    
This is (will be) fixed in PostGIS 2.0.1: trac.osgeo.org/postgis/ticket/1782 –  Steve Pike Jun 14 '12 at 10:34

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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