How do I get the area of a GeoJSON polygon with Python

Question

I have the GeoJSON

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [[13.65374516425911, 52.38533382814119], [13.65239769133293, 52.38675829106993], [13.64970274383571, 52.38675829106993], [13.64835527090953, 52.38533382814119], [13.64970274383571, 52.38390931824483], [13.65239769133293, 52.38390931824483], [13.65374516425911, 52.38533382814119]]
        ]
      }
    }
  ]
}

which http://geojson.io displays as

I would like to calculate its area (87106.33m^2) with Python. How do I do that?

What I tried

# core modules
from functools import partial

# 3rd pary modules
from shapely.geometry import Polygon
from shapely.ops import transform
import pyproj

l = [[13.65374516425911, 52.38533382814119, 0.0], [13.65239769133293, 52.38675829106993, 0.0], [13.64970274383571, 52.38675829106993, 0.0], [13.64835527090953, 52.38533382814119, 0.0], [13.64970274383571, 52.38390931824483, 0.0], [13.65239769133293, 52.38390931824483, 0.0], [13.65374516425911, 52.38533382814119, 0.0]]
polygon = Polygon(l)

print(polygon.area)
proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),
                   pyproj.Proj(init='epsg:3857'))
print(transform(proj, polygon).area)

It gives 1.1516745933889345e-05 and 233827.03300877335 - that the first one doesn't make any sense was expected, but how do I fix the second one? (I have no idea how to set the pyproj.Proj init parameter)

I guess epsg:4326 makes sense at it is WGS84 (source), but for epsg:3857 I'm uncertain.

Better results

The following is a lot closer:

# core modules
from functools import partial

# 3rd pary modules
import pyproj
from shapely.geometry import Polygon
import shapely.ops as ops

l = [[13.65374516425911, 52.38533382814119, 0],
     [13.65239769133293, 52.38675829106993, 0],
     [13.64970274383571, 52.38675829106993, 0],
     [13.64835527090953, 52.38533382814119, 0],
     [13.64970274383571, 52.38390931824483, 0],
     [13.65239769133293, 52.38390931824483, 0],
     [13.65374516425911, 52.38533382814119, 0]]
polygon = Polygon(l)

print(polygon.area)
geom_area = ops.transform(
    partial(
        pyproj.transform,
        pyproj.Proj(init='EPSG:4326'),
        pyproj.Proj(
            proj='aea',
            lat1=polygon.bounds[1],
            lat2=polygon.bounds[3])),
    polygon)
print(geom_area.area)

it gives 87254.7m^2 - that is still 148m^2 different from what geojson.io says. Why is that the case?

Oculist

Because Web Mercator doesn't preserve area: geography.hunter.cuny.edu/~jochen/GTECH361/lectures/lecture04/…. Welcome to map projections. Your question is a better for Geographic Information Systems (There's nothing wrong with your code, just your choice of projection.) and is better off being specific about getting the wrong value for the area. And here's a link to the actual 3857 projection; the one you link is something else. – Aggravation

For anyone trying our the attached python code and getting an error, it should be lat_1 and lat_2 instead of lat1 and lat2 – Thickset

Answer 1

It looks like geojson.io is not calculating the area after projecting the spherical coordinates onto a plane like you are, but rather using a specific algorithm for calculating the area of a polygon on the surface of a sphere, directly from the WGS84 coordinates. If you want to recreate it you can find the source code here.

If you are happy to carry on projecting the coordinates to a flat system to calculate the area, since it's good enough accuracy for your use case, then you might trying using this projection for Germany instead. E.g:

from osgeo import ogr
from osgeo import osr

source = osr.SpatialReference()
source.ImportFromEPSG(4326)

target = osr.SpatialReference()
target.ImportFromEPSG(5243)

transform = osr.CoordinateTransformation(source, target)

poly = ogr.CreateGeometryFromJson(str(geoJSON['features'][0]['geometry']))
poly.Transform(transform)
poly.GetArea()

which returns 87127.2534625642

Telex

Hmm... Is there a simple answer what the correct projection is? What are the main differences between 4326 and 5243? – Oculist

@MartinThoma There is no such thing as a simple answer to "the correct projection." All the projections have different properties. It very much depends on your specific use case. In particular, a lot of projections give accurate values only in certain areas of the world. – Aggravation

This is super interesting! Could you name two reasonable projections for Berlin and their advantages compared to each other? – Oculist

Answer 2

There is a port of Mapbox's geojson-area for Python: https://pypi.org/project/area/

from area import area

geoJSON = json.loads(...)

print(area(geoJSON['features'][0]['geometry']))

With your GeoJSON, it returns 87106.33 like on http://geojson.io/.

Indirection

Answer 3

If using GeoDjango, as explained here, you can project the geojson first and retrieve the area in square meters from that:

def get_acres(self): 
    """ 
    Returns the area in acres. 
    """ 
    # Convert our geographic polygons (in WGS84)
    # into a local projection - for New York here as an example (EPSG:32118)
    self.polygon.transform(32118)
    meters_sq = self.polygon.area
    
    acres = meters_sq * 0.000247105381  # meters^2 to acres

    return acres

Sherrie