0

I have a question with a query in PostGIS. I try to calculate the area from a polygon. It actually works but the area-result is not correct. My Query is:

SELECT
name, ST_Area(way)/10000 AS hectares
FROM
osm_data
ORDER BY hectares DESC
LIMIT 1;

The problem is, that the amount of ha is completely wrong. So i tryed to change the EPSG in the standard 4326 EPSG. But when I use than the query it shows me an error:

SELECT
name,ST_Transform ((way, 4326)ST_Area(way))/10000 AS hectares
FROM
planet_osm_polygon
ORDER BY hectares DESC
LIMIT 1;

SQL state: 42601

Where is the mistake?!

Vince
  • 20,017
  • 15
  • 45
  • 64
Ayron
  • 43
  • 2
  • 7
  • 2
    If the data has an SRID or 4326, your area is calculated in Cartesian degrees. Your second attempt has a clear syntax error. You can cast to geography or provide a proper syntax for the ST_Transform (placing it inside the ST_Area) – Vince Sep 27 '18 at 15:42

2 Answers2

6

OSM data is in webmercator projection by default (EPSG 3857), which is a good projection for display but that introduces severe distortions (about 40% at 45 degrees of latitude). You can look at this article for the explanations.

To compute the distance, you have two choices:

1) you can transform the data to a projection that is suitable to your area and that does not distort the data too much (UTM or alike might be good candidates).

2) you transform your data to geography, meaning it is not projected on a flat ground, but it is rather kept on a round earth surface. This approach is easy but it considers great-arcs instead of straight lines between two points. Unlike 4326, the unit of the geography type is meters, not degrees.

The doc says

The basis for the PostGIS geometry type is a plane. The shortest path between two points on the plane is a straight line. That means calculations on geometries (areas, distances, lengths, intersections, etc) can be calculated using cartesian mathematics and straight line vectors.

The basis for the PostGIS geographic type is a sphere. The shortest path between two points on the sphere is a great circle arc. That means that calculations on geographies (areas, distances, lengths, intersections, etc) must be calculated on the sphere, using more complicated mathematics. For more accurate measurements, the calculations must take the actual spheroidal shape of the world into account, and the mathematics becomes very complicated indeed.

Assuming your data is in EPSG 3857, you could compute the area using:

SELECT
name, ST_Area(ST_Transform(way, 4326)::geography) AS m2,  
      ST_Area(ST_Transform(way, 4326)::geography)/10000 as ha
FROM
planet_osm_polygon
ORDER BY ha DESC
LIMIT 1;
JGH
  • 41,794
  • 3
  • 43
  • 89
  • I didn't know that that was possible to cast output to spherical distance units with the ::geography suffix. Like I mentioned in my answer, my PostGIS is a bit rusty. Interesting. – Trevor J. Smith Sep 27 '18 at 18:27
  • @TrevorJ.Smith it works only if the source SRID is unknown (not set) or 4326 – JGH Sep 27 '18 at 18:32
  • @Trevor it assumes EPSG:4326 even if not set, so having SRID:0 is a possible source of error. – geozelot Sep 27 '18 at 19:27
  • Thanks, for the explanation and the query both makes completely sense! Query works as well! – Ayron Oct 02 '18 at 08:35
-1

I may be a bit rusty, but I think your problems are that the syntax is wrong (it looks like you're calling area as a method, which isn't a thing in SQL) and that you aren't nesting the transform function appropriately. You need to perform the transform BEFORE you calculate your area:

SELECT
name, ST_Area(ST_Transform(way, 4326)) AS squared_degrees
FROM
planet_osm_polygon
ORDER BY hectares DESC
LIMIT 1;

The ST_Area returns a finite number, not a geometry, and so it isn't possible to perform a spatial transform on it.

Seeing as you're a new user, please let us know if this helps by accepting my answer, updating your question with more information or commenting on my answer if this seems unclear. Best of luck!

EDIT: In order to convert this unit to hectares, you need to cast this spatially transformed unit to a geography with arc-lengths (ST_Transform(way, 4326)::geography) and then divide that number by 10000. See JGH's answer for more information.

Trevor J. Smith
  • 1,398
  • 11
  • 26