CODEWITHSUNDEEP
pythongeojsonpyproj
Martin Thoma
Martin Thoma

Reputation: 136237

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

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

enter image description here

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?

Upvotes: 13

Views: 20679

Answers (3)

Akaisteph7
Akaisteph7

Reputation: 6476

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

Upvotes: 0

Benedikt Köppel
Benedikt Köppel

Reputation: 5109

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/.

Upvotes: 5

fordy
fordy

Reputation: 2691

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

Upvotes: 11

Related Questions