Ben
Ben

Reputation: 71

I tried all ways, but still my area is calculated wrongly in Postgis

I created a very simple polygon in the middle of Germany to demonstrate my problem.

You can visualize it in geojsonlint using the following GeoJSON

{"type":"Polygon","coordinates":[[
[10.439844131469727,51.17460781257472],
[10.430574417114258,51.1753073564544],
[10.429565906524658,51.17179607723465],
[10.438792705535889,51.170706315523866],
[10.439372062683105,51.17267055874809],
[10.43975830078125,51.17439256616884],
[10.439844131469727,51.17460781257472]]]G}

When calculating the surface with online tools (e.g. http://www.daftlogic.com/projects-google-maps-area-calculator-tool.htm, but I tried several), I get the following numbers (these are based on a similar drawing of the polygon, but not the exact same one, as I couldn't copy it over to these tools):

Now I want to calculate these areas using POSTGIS, but I always get wrong and not matching numbers.

First I started without transformation using the examples given here:

http://postgis.net/docs/ST_Area.html

 SELECT ST_Area(the_geom) As sqft, ST_Area(the_geom)*POWER(0.3048,2) As sqm
 FROM (SELECT ST_GeomFromText('
 POLYGON ((51.17460781257472  10.439844131469727,
 51.1753073564544 10.430574417114258, 
 51.17179607723465 10.429565906524658,
 51.170706315523866 10.438792705535889,
 51.17267055874809 0.439372062683105,
 51.17439256616884 10.43975830078125,
 51.17460781257472 10.439844131469727))',4326) ) As foo(the_geom);

--> sqft = 3.52643124351653e-05 and sqm = 3.27616182873666e-06

How can I interprete these numbers? Then I tried to transform it to WGS 84 / UTM zone 33N 32633

SELECT ST_Area(the_geom) As sqft, ST_Area(the_geom)*POWER(0.3048,2) As sqm
FROM (SELECT ST_Transform(ST_GeomFromText('
POLYGON ((51.174661624019286 10.440187454223633,
51.17067940750161 10.438899993896484,
51.17197097486416 10.429544448852539, 
51.17536116708255 10.430488586425781,
51.174661624019286 10.440187454223633))',4326),32633) ) As foo(the_geom);

--> sqft = 662918.939349234 and sqm = 61587.1847391195

But even these numbers don't come close.

Upvotes: 0

Views: 2944

Answers (2)

John Powell
John Powell

Reputation: 12581

I converted the coordinates into EPSG: 31467, see epsg:31467 which is projected to meters and applies to the area of Germany covered by your geometry.

select st_area(st_transform(st_setsrid(st_geomfromtext('POLYGON((10.439844131469727    
  51.17460781257472,10.430574417114258 51.1753073564544,10.429565906524658
  51.17179607723465,10.438792705535889 51.170706315523866, 10.439372062683105 
  51.17267055874809, 10.43975830078125 51.17439256616884, 10.439844131469727 
  51.17460781257472))'),4326),31467));

and got the answer: 274442.27 m2 which is within 0.007% of your original answer.

Measurements are usually more accurate in projected coordinate systems that use a geoid appropriate to that region. If you run this query on the spatial reference system table in Postgis for that projection:

 select * from spatial_ref_sys where srid=31467;

you will see some more details, such as the fact that it uses the Bessel 1841 spheroid.

EDIT: your original geojson has coordinates in x/y, but for some reason you flipped them when putting them into Postgis.

Upvotes: 1

Ben
Ben

Reputation: 71

The coordinates of the polygon were accidentally loaded as lat,lon instead of lon, lat.

http://postgis.net/2013/08/18/tip_lon_lat says

In spatial databases spatial coordinates are in x = longitude, and y = latitude

Upvotes: 2

Related Questions