Opened 15 years ago
Closed 15 years ago
#304 closed defect (fixed)
Geography Best SRID or ?
Reported by: | robe | Owned by: | pramsey |
---|---|---|---|
Priority: | medium | Milestone: | PostGIS 1.5.0 |
Component: | postgis | Version: | master |
Keywords: | Cc: |
Description
I'm a bit puzzled maybe its because I haven't thought about how the curvature of the earth affects area much, but I would expect that if I buffer a point, any point on the globe by say 10 meters, that each polygon would be the same area. It more or less is if I convert to UTM zone, buffer, convert back to geography —> take the area. But not when I use ST_Buffer geography. What gives?
(If I buffer by 10 meters — I expect the answer to be pi()*102 (314 +- 2). Really 312 because of the inaccuracy of the default buffer (using 8 segs per quarter circle).
and it is for the most part except when I go to Helsinki, Bergen. My plot of land is around 288 sq meters. I guess I shouldn't be visiting those places.
SELECT city, ST_Area(ST_Buffer(geography(the_geom), 10),true) As geog_area, ST_Area(geography(ST_Transform(ST_Buffer(ST_Transform(the_geom,utm_srid), 10), 4326))) As utm_geog_area FROM (VALUES ('Helsinki', ST_GeomFromText('POINT(24.9767 60.1964)',4326), 32635 ), ('Bergen', ST_GeomFromText('POINT(5.4907 60.35)',4326), 32631 ) ) f(city, the_geom,utm_srid);
the above query yields
city | geog_area | utm_geog_area ---------+------------------+------------------ Helsinki | 287.711725821719 | 312.298024061475 Bergen | 288.120740812147 | 312.249761770185
Change History (5)
comment:1 by , 15 years ago
Resolution: | → fixed |
---|---|
Status: | new → closed |
comment:2 by , 15 years ago
Resolution: | fixed |
---|---|
Status: | closed → reopened |
Not so fast cowboy. cough cough. While it now passes the 2 above tests, it still fails my more exhaustive shotgun test even in cases where a simple utm hack gives a pleasing answer.
--helper function copied from wiki CREATE OR REPLACE FUNCTION utmzone(geometry) RETURNS integer AS $BODY$ DECLARE geomgeog geometry; zone int; pref int; BEGIN geomgeog:= ST_Transform($1,4326); IF (ST_Y(geomgeog))>0 THEN pref:=32600; ELSE pref:=32700; END IF; zone:=floor((ST_X(geomgeog)+180)/6)+1; RETURN zone+pref; END; $BODY$ LANGUAGE 'plpgsql' IMMUTABLE COST 100; -- shotgun test SELECT ST_AsText(the_geog) as the_pt, ST_Area(ST_Buffer(the_geog,10)) As the_area, ST_Area(geography(ST_Transform(ST_Buffer(ST_Transform(geometry(the_geog),utm_srid),10),4326))) As geog_utm_area FROM (SELECT geography(ST_SetSRID(ST_Point(i*0.5,j*0.5),4326)) As the_geog, utmzone(ST_SetSRID(ST_Point(i*0.5,j*0.5),4326)) As utm_srid FROM generate_series(-350,350) As i CROSS JOIN generate_series(-175,175) As j ) As foo WHERE ST_Area(ST_Buffer(the_geog,10)) NOT between 310 and 314 LIMIT 10;
Do not pass go until you can pass the above (e.g. no answers returned), or utm doesn't give a better answer or you can give a 4 paragraph or less reason why your best is the best you can do.
comment:3 by , 15 years ago
Well, there was a boneheaded mistake that was using Arctic stereographic in Antarctica. Using the correct stereographic fixes that, at r4822. I still don't pass your test though, perhaps because stereographic isn't All That? The trouble w/ UTM at the poles is that the "right" zone to use gets very indeterminate.
comment:4 by , 15 years ago
Latitude of true scale for Polar Stereographic is at 70S http://www.classroom.antarctica.gov.au/introduction/references-and-resources/classroom-antarctica-resources/maps/map-projections-in-the-antarctic-and-their-use/ which is coincidentally where the calculations agree.
Those northerly cities are falling into the polar stereographic zone, which perhaps runs too far south. I've moved the cutoff line another five degrees north to 65. r4812