11

I need to calculate areas and intersection areas for polygons (some real geographical objects like lake, city, country, e.t.c.). Polygons located in California, New Zealand, Russia.Anadyr, Sweden

All polygons are in WGS84.

What I did using GeoTool java API:

  1. Project all polygons using EPSG:3488, EPSG:NAD83(NSRS2007) / California Albers and calculated areas and overlap areas.
  2. Did the same using World_Mollweide and World_Eckert_IV
  3. Picked "local specific projections" for polygons from California, New Zealand, e.t.c.

I assume that #3 is the most accurate result, since I pick projection which covers polygon area

Result:

'#2 showed the worst result comparing to #3

'#1 and #3 areas and intersection areas difference is less than 0.1%

Why? I pick absolutely wrong projection EPSG:3488 (California) for polygons from Sweden and get pretty the same areas and intersection areas?

UPD: Looks like I didn't explain my confusion correctly. Here is sample output with explanation

#area_from_new_zealand_1
EPSG_27200 area[11733479] CRS[World_Mollweide] area[11736023] diff[2544] [0.0%]
EPSG_27200 area[11733479] CRS[World_Eckert_IV] area[11736033] diff[2554] [0.0%] 
EPSG_27200 area[11733479] CRS[EPSG:NAD83(NSRS2007) / California Albers] area[11736034] diff[2555] [0.0%]

#area_from_new_zealand_2
EPSG_27200 area[2952725]  CRS[World_Mollweide] area[2953281] diff[556] [0.0%] 
EPSG_27200 area[2952725]  CRS[World_Eckert_IV] area[2953342] diff[617] [0.0%] 
EPSG_27200 area[2952725]  CRS[EPSG:NAD83(NSRS2007) / California Albers] area[2953467] diff[743] [0.0%]


#intersection_area_between_two_new_zealand_areas
EPSG_27200 intersection area[1001857] CRS[World_Mollweide]                          area[1002082] diff[225] [0.0%] 
EPSG_27200 intersection area[1001857] CRS[World_Eckert_IV]                          area[1002082] diff[225] [0.0%] 
EPSG_27200 intersection area[1001857] CRS[EPSG:NAD83(NSRS2007) / California Albers] area[1002096] diff[239] [0.0%]


#area_from_alaska_1
EPSG_3338 area[56278347]    CRS[World_Mollweide] area[56041510] diff[236837] [0.4%] 
EPSG_3338 area[56278347]    CRS[World_Eckert_IV] area[56041585] diff[236763] [0.4%] 
EPSG_3338 area[56278347]    CRS[EPSG:NAD83(NSRS2007) / California Albers] area[56278426] diff[79] [0.0%]


#area_from_alaska_2
EPSG_3338 area[17564799282] CRS[World_Mollweide] area[17486015889] diff[78783393] [0.4%] 
EPSG_3338 area[17564799282] CRS[World_Eckert_IV] area[17486869816] diff[77929466] [0.4%]
EPSG_3338 area[17564799282] CRS[EPSG:NAD83(NSRS2007) / California Albers] area[17566197286] diff[1398004] [0.0%]


#intersection_area_between_two_alaska_areas 
EPSG_3338 intersection area[43808167] CRS[World_Mollweide] area[45066901] diff[1258734] [2.8%] 
EPSG_3338 intersection area[43808167] CRS[World_Eckert_IV] area[45163183] diff[1355016] [3.0%] 
EPSG_3338 intersection area[43808167] CRS[EPSG:NAD83(NSRS2007) / California Albers] area[43885182] diff[77015] [0.2%]

My confusion is: EPSG:3488 designed to be used in California

I pick "wrong" projection EPSG:3488 for Alaska, New Zealand areas and see that resulting calculations don't differ "significantly" from correct projections. EPSG:3488 even performs better than Mollweide, Eckert_IV projections designed to be used around the world.

Community
  • 1
Capacytron
  • 327
  • 2
  • 11
  • I have also found that there is close to no observable difference between these two projections, however the difference still exists. In ArcGIS, you cannot create a "feature dataset" unless your data is in the same projection even with such a small difference as found between WGS84 and NAD83. The following webpage was very informative to me and I hope you find it to be useful as well. http://www.differencebetween.net/technology/difference-between-wgs84-and-nad83/ I would have put this as a comment but I do not have 50 rep :( – Justin Q Aug 13 '15 at 16:02
  • what are you comparing the results to? – Ian Turton Aug 14 '15 at 09:07
  • @iant please see updated question. I added comparison output. – Capacytron Aug 14 '15 at 09:37
  • You could try the AUTO projections (UTM centred on a user supplied point) - construct CRS as String code = "AUTO:42001," + x + "," + y; // System.out.println(code); CoordinateReferenceSystem auto = CRS.decode(code); – Ian Turton Aug 14 '15 at 09:50

2 Answers2

9

"EPSG:3488, EPSG:NAD83(NSRS2007) / California Albers" is an equal-area projection. It is based on the Albers Conic, which is defined for the northern hemisphere. Because Sweden is within its range of definition, it is equal-area in Sweden. This means that (up to floating point rounding error) it will give absolutely correct areas.

Neither the Mollweide nor the Eckert are exactly equal-area, but (as M. Kennedy kindly points out in a comment) they are approximately so. The distortions they introduce will be comparable to the differences between the sphere and the ellipsoid, which are limited to about one part in 300 (0.3%).

Community
  • 1
whuber
  • 69,783
  • 15
  • 186
  • 281
  • 1
    Bill, Eckert IV and Mollweide are equal-area projections, but only have spherical algorithms. – mkennedy Aug 13 '15 at 16:41
  • 1
    @mkennedy Oops--I should have checked. Thank you for the correction, Melita. I'll fix that comment. – whuber Aug 13 '15 at 16:41
  • 1
    @mkennedy so ESRI:53009 and ESRI:54009 are actually identical? I see that Snyder does not give formulas on the ellispoid, so what's the sense of those World_xxx projections from ESRI based on WGS84? – AndreJ Aug 13 '15 at 17:14
  • 1
    Esri:53009 (and the other 53xxx entries) use a sphere-based GeoCRS with R=6371000.0 m. The Esri:54xxx range uses a WGS84 geoCRS so the actual radius used is the semimajor axis, 6378137.0. They were both added just as test/sample definitions. – mkennedy Aug 13 '15 at 17:25
  • 1
    @Andrej For area calculations, one can always compute the authalic latitudes and apply the projection formulas to them if that kind of accuracy is needed. – whuber Aug 13 '15 at 19:31

  • it is equal-area in Sweden. This means that (up to floating point rounding error) EPSG:3488, shows 0.0% diff. It's more close to south pole than to north. Sorry it it's stupid question, I'm developer, have no ph.d in Geography.

     – Capacytron Aug 13 '15 at 22:06
  • 2
    The South Pole is tens of thousands of kilometers from Sweden: it's in Antarctica. The North pole is only a few thousand kilometers from Sweden. But it doesn't matter: if you use an equal area projection and it is capable of projecting a region whose area you want to compute, then it will compute the area correctly (for the datum on which it is based). Some equal-area projections are capable of projecting the entire world except for a single point (which depends on the projection). – whuber Aug 13 '15 at 22:12
  • @whuber, one of my comments didn't post correct. I tested New Zealand areas with EPSG:3488 and EPSG_27200. They show < 0.0% diff in area calculation and area intersection calculation. New Zealand is close to South Pole and I still get correct result. – Capacytron Aug 14 '15 at 07:30
  • You will get exact results with any equal area projections except for one point, and the worldwide projections that use a sphere (although World_Mollweide claims to use the WGS84 ellipsoid, but as mkennedy stated that is not true) – AndreJ Aug 14 '15 at 09:52
  • @AndreJ (This is just nit-picking, but the idea might be of interest to readers.) You're correct about the "any" part, but not about the "except for one point": not all equal-area projections are defined that extensively. For instance, it's mathematically impossible for any cylindrical equal-area projection to be defined at its two poles. Oblique aspects of most worldwide ellipsoidal projections (of any sort) will also fail along a strange-looking curve near the diametrically opposite point of their origins. – whuber Aug 14 '15 at 13:36
1

@whuber's assertion that an equal-area projection "will give absolutely correct areas" comes with an asterisk, namely, assuming that the edges of the polygon are straight lines in said projection. This is often a good approximation, particularly if the edges are short; but it is rarely strictly true.

If, on the other hand, the edges of your polygon are geodesics or rhumb lines, other techniques can be used to determine the area accurate to round-off. My online planimeter implements these. Give it a try.

cffk
  • 3,271
  • 18
  • 24
  • Hi! Thanks for the all input. So want could be the summary? EPSG:3488 uses Albers equal-area projection, that is why it correctly calculates areas all around the world, even on South pole? – Capacytron Aug 18 '15 at 08:54
  • Probably Albers equal-area is going to provide a decent enough result for your application. However blindly using this projection to calculate the area of Antarctica (which encircles the pole) will give a nonsense result. So for general use, I would recommend calculating the geodesic area so you don't have to worry about the limitations of each particular equal-area projection. – cffk Aug 18 '15 at 09:49
  • Thanks for the reply. Unfortunately, we need square area in meters. – Capacytron Aug 18 '15 at 09:56
  • The Planimeter tool gives you the area in square meters. – cffk Aug 18 '15 at 10:58
  • The same functionality is available as a Java package, Geographiclib-Java. The documentation includes code to compute the area of a polygon specified as a set of latitude/longitude points with the result given in square meters. – cffk Aug 18 '15 at 11:09
  • we are using geotools to project WGS to planar and then cal area and intersections using JTS. You suppose that Geographiclib-Java could give better performance? – Capacytron Aug 18 '15 at 11:13
  • GeographicLib-Java obviates the need to project the data to a planar system and it calculates accurately the area of a polygon whose edges are geodesics. It does not deal with intersections of polygons. – cffk Aug 18 '15 at 11:29