Shapely: Polygon area calculation doesn't match other tools like google earth engine or geojson.io

[brentonmallen1]

Polygon area calculation doesn't seem to match up with other similar tools such as Geojson.io or Google Earth Engine. For example:

coords = [(-97.59238135821987, 43.47456565304017),
 (-97.59244690469288, 43.47962399877412),
 (-97.59191951546768, 43.47962728271748),
 (-97.59185396090983, 43.47456565304017),
 (-97.59238135821987, 43.47456565304017)]

projection = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'), pyproj.Proj(init='epsg:3857'))
shapely.ops.transform(projection, shapely.geometry.Polygon(sample_coords)).area

yields an area of 45573.993884405005 m^2 while Google Earth Engine and Geojson.io report 23944.14737277293 and 23997.77, respectively.

The area in shapely is calculated here quoting a method whose link is broken, while the Geojson.io is calculated here via a method described on starting on page 3 of this document.

I've converted this briefly to python with the argument expecting the same coordinate scheme as above (list of coordinates):

def ringArea(coords):
    #var p1, p2, p3, lowerIndex, middleIndex, upperIndex, i,
    area = 0
    wgs84_radius = 6378137 # earth's equitorial radius in wgs84 = espg 4326
    coordsLength = len(coords)

    if coordsLength > 2:
        for i in range(coordsLength):
            if i == coordsLength - 2:
                lowerIndex = coordsLength - 2
                middleIndex = coordsLength - 1
                upperIndex = 0
            elif i == coordsLength - 1:
                lowerIndex = coordsLength - 1
                middleIndex = 0
                upperIndex = 1
            else:
                lowerIndex = i
                middleIndex = i + 1
                upperIndex = i + 2

            p1 = coords[lowerIndex]
            p2 = coords[middleIndex]
            p3 = coords[upperIndex]
            area += (np.deg2rad(p3[0]) - np.deg2rad(p1[0]) ) * np.sin(np.deg2rad(p2[1]))

        area = area * wgs84_radius * wgs84_radius / 2
    return area

This method returns an area of 23997.769250450423 m^2.

My main question is: Is the current method preferred for some reason? If so, can the link to the cited be fixed (i'm not sure where it's supposed to go). If not, is it possible to conform (or add the option) to the method described int he JPL doc?

I'd be happy to contribute to this if the option for this other method would like to be added but I didn't see any contributing guidelines.

[dnomadb]

As far as I understand, shapely assumes all points are cartesian, and calculates area in whatever coordinates it is given; in the example above wgs coordinates are converted to web mercator epsg3857 to calculate area, but web mercator is not an equal area projection.

If we substitute conversion to an equal area coordinate system such as world mollweide (https://spatialreference.org/ref/esri/54009/) , the area calculation is performed as expected:

coords = [(-97.59238135821987, 43.47456565304017),
 (-97.59244690469288, 43.47962399877412),
 (-97.59191951546768, 43.47962728271748),
 (-97.59185396090983, 43.47456565304017),
 (-97.59238135821987, 43.47456565304017)]

projection = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'), pyproj.Proj(init='esri:54009'))
shapely.ops.transform(projection, shapely.geometry.Polygon(coords)).area

Which gives 23997.775332830173.

We use this over in https://github.com/mapbox/fio-area. Hope this helps!

[mwtoews]

Shapely only does Cartesian calculations, so this issue is mostly off topic.

I've just fixed the broken link in shapely/algorithms/cga.py to https://web.archive.org/web/20080209143651/http://cgafaq.info:80/wiki/Polygon_Area

There are several different ways to calculate areas from geodetic coordinates, and they will each give a slightly different result. To add another method, use Planimeter on WGS84 for the example to get 23988.9 m^2 with an error about 0.2 m^2, so I'd call this the "best answer". This method is also used in PostGIS on geography types.

[sgillies]

Thanks @dnomadb @mwtoews !

[brentonmallen1]

I appreciate the input from everyone. I'm pretty new to the GIS world so it's all helpful, thank you!

[brentonmallen1]

@dnomadb when trying the method suggested I got the following error:

CRSError: Invalid projection: +init=epsg:54009 +type=crs: (Internal Proj Error: proj_create: crs not found)

Have you come across this?

[dnomadb]

@brentonmallen1 it looks like you have epsg:54009 where it should be esri:54009 (as in example here: https://github.com/Toblerity/Shapely/issues/726#issuecomment-499252744 -- it's an easy thing to mix-up)

[brentonmallen1]

That was a typo on my part trying different things to see if I could get it to work. I still get the error with esri:54009:

projection = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'), pyproj.Proj(init='esri:54009'))

```---------------------------------------------------------------------------
CRSError Traceback (most recent call last)
in
5 (-97.59238135821987, 43.47456565304017)]
6
----> 7 projection = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'), pyproj.Proj(init='esri:54009'))
8 shapely.ops.transform(projection, shapely.geometry.Polygon(coords)).area

~/miniconda3/envs/fs-ml/lib/python3.6/site-packages/pyproj/proj.py in __init__(self, projparams, preserve_units, **kwargs)
145 '116.366 39.867'
146 """
--> 147 self.crs = CRS.from_user_input(projparams if projparams is not None else kwargs)
148 # make sure units are meters if preserve_units is False.
149 if not preserve_units and "foot" in self.crs.axis_info[0].unit_name:

~/miniconda3/envs/fs-ml/lib/python3.6/site-packages/pyproj/crs.py in from_user_input(cls, value)
389 if isinstance(value, _CRS):
390 return value
--> 391 return cls(value)
392
393 def get_geod(self):

~/miniconda3/envs/fs-ml/lib/python3.6/site-packages/pyproj/crs.py in __init__(self, projparams, **kwargs)
258 raise CRSError("Invalid CRS input: {!r}".format(projparams))
259
--> 260 super(CRS, self).__init__(projstring)
261
262 @classmethod

pyproj/_crs.pyx in pyproj._crs._CRS.__init__()

CRSError: Invalid projection: +init=esri:54009 +type=crs: (Internal Proj Error: proj_create: cannot expand +init=esri:54009 +type=crs)```

[brentonmallen1]

In case anyone else comes across the above issue w.r.t. re-projecting to ESRI:54009, here is a solution https://github.com/pyproj4/pyproj/issues/341