Geography Best SRID or ?

[robe]

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()*10^{2 (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

[pramsey]

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

[robe]

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.

[pramsey]

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.

[pramsey]

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.

[pramsey]

King me. r4823