How do I get the area of a GeoJSON polygon with Python
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
Martin Thoma
I also have a blog about Code, the Web and Cyberculture (medium as well) and a career profile on Stackoverflow. My interests are mainly machine-learning, neural-networks, data-analysis, python, and in general backend development. I love building secure and maintainable systems.
Updated on June 04, 2022Comments
-
Martin Thoma almost 2 years
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
and233827.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 thepyproj.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?