0

I'm working in the SRID 4326 system and need to compute the area of rectilinear polygons. Our geometry (geom) column is in WKB format

When I execute the ST_Area function it appears my results are returned in squared degrees when I execute a call as follows:

public.ST_Area(geom)/1000000 from production.schemaA where item = '123456'

Then if I cast the geom column to geography the result returned was 33.97 km2

This is still quite a bit off from what I get checking the area of the SAME polygon in Arc Map and QGIS Desktop (Field calculator method) about 43.37 km2

I have investigated using the ST_Transform function to transform out of 4326 but this doesn't work as you have to know which UTM/WGS SRID to use in advance based on the (UTM/ZONE) zone you are working in. We'll be hitting our global database of polygons all over the world had identifying the specific zone is cumbersome.

Example: public.ST_Area(geom::geography)/1000000 from production.schemaA where item = '123456'

There are threads that are close to this issue but not exactly the same issue of have no solution.

  • 3
    Cast to geography should give correct area. Could you add one geometry as WKT and areas which you get from PostGIS, QGIS, and ArcMap? – user30184 Aug 31 '18 at 12:26
  • 'POLYGON ((38.31000000 36.85999999, 38,31000000 36.90999999, 38,38000000 36.90999999, 38.38000000 36.8599999)}' – user9491577 Aug 31 '18 at 12:38
  • That's not a valid polygon. Check the decimal separators and make sure that first and last points are the same. – user30184 Aug 31 '18 at 12:46
  • My bad -- typed comma as instead of dot: 'POLYGON ((38.31000000 36.85999999, 38.31000000 36.90999999, 38.38000000 36.90999999, 38.38000000 36.85999999))' – user9491577 Aug 31 '18 at 12:51
  • It was still not closed but I think that you mean this POLYGON (( 38.31 36.85999999, 38.31 36.90999999, 38.38 36.90999999, 38.38 36.85999999, 38.31 36.85999999 )) – user30184 Aug 31 '18 at 12:57
1

Compare the results from PostGIS with geography:

select ST_Area(
ST_GeogFromText(
'POLYGON ((38.31000000 36.85999999, 38.31000000 36.90999999, 38.38000000 36.90999999, 38.38000000 36.85999999, 38.31000000 36.85999999))'));

34625394.4708328

and with geometry transformed into UTM zone 37 (EPSG:32637):

select ST_Area(
ST_Transform(
ST_GeomFromText(
'POLYGON ((38.31000000 36.85999999, 38.31000000 36.90999999, 38.38000000 36.90999999, 38.38000000 36.85999999, 38.31000000 36.85999999))',4326),32637));


34600605.4910705

Results are about the same, taking into account that UTM has a scale factor of 0.9996. It is unclear in which projection QGIS and ArcMap are computing the areas but the results do not seem to fit well with the ground truth. That's a drawback of software which do clever on-the-fly transformations which can make it hard for the users to know what really happens.

  • The problem with method one is you have to know the zone prior and this would have to be determined by the application before the DB select statement is invoked. That over complicates the task of retrieving the the area on the backend. So for our implementation the second method is the way to go. I'm a bit concerned about the area difference I get via QGIS Desktop and Arcmap versus the areas retrieved above though, – user9491577 Aug 31 '18 at 14:10
  • It's a bit of a hack, but, assuming your coordinates will always fit into one UTM zone, you could always substract the minimum y coordinate from all the coordinates, forcing the polygon into UTM zone 31, ie, with London on the Western edge, and then use the area for that SRID, with the same 0.9996 scale factor proviso as before. The area will not change by changing the longitude by a constant amount, after all. – John Powell Aug 31 '18 at 14:16
  • Getting back to this query. I trying to execute the following SQL command based on your input::SELECT stockId,public.ST_Area( public.ST_GeogFromText(geom)) FROM temporary.myschema. I'm getting a error attempting to resolve the the geom column which in our case is polygons in WKB format. Is the issue related to our geom column is WKB and I'm using ST_GeogFromText. I also tried ST_GeomFromText with no success. – user9491577 Sep 6 '18 at 11:55
  • I modified the above query as follows: SELECT stockId,public.ST_Area( public.ST_GeogFromText(public.ST_AsText(geom))) FROM temporary.myschema and achieved results. – user9491577 Sep 6 '18 at 12:03
  • I believe that whats i was missing -- the conversion to WKT – user9491577 Sep 6 '18 at 12:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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