13

I have polygons from coordinates in (python shapely) that looks like this

POLYGON ((24.8085317 46.8512821, 24.7986952 46.8574619, 24.8088238 46.8664741, 24.8155239 46.8576335, 24.8085317 46.8512821))

I would like to calculate the area of this polygon in km^2. What would be the best way to do this in Python?

  • 1
    You may look at stackoverflow.com/questions/23697374/… – SIslam Dec 26 '14 at 9:12
  • If you're getting the following error implementing one of the above solutions, is because the lat1 and lat2 should be lat_1 and lat_2: pyproj.exceptions.CRSError: Invalid projection: +proj=aea +lat1=37.843975868971484 +lat2=37.844325658890924 +type=crs: (Internal Proj Error: proj_create: Error -21: conic lat_1 = -lat_2) – Ramtin Kermani Apr 17 at 0:18
15

It looks like your coordinates are longitude and latitude, yes? Use Shapely's shapely.ops.transform function to transform the polygon to projected equal area coordinates and then take the area.

python
import pyproj
from functools import partial

geom_aea = transform(
partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj(
        proj='aea',
        lat1=geom.bounds[1],
        lat2=geom.bounds[3])),
geom)

print(geom_aea.area)
# Output in m^2: 1083461.9234313113 
  • You should probably indicate that partial is not a built-in; pyproj will have to be imported and possibly installed, etc. – kingledion Apr 26 '18 at 18:27
  • I noticed pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2) resulted in CRSError: Invalid projection: +proj=aea +lat1=5.0 +lat2=6.0 +type=crs. Changing lat{1,2} into lat_{1,2} as implied by proj4 documentation fixed it: pyproj.Proj(proj='aea', lat1=lat1, lat2=lat2). Is this Answer accurate, or should it be updated? – Herbert Mar 18 at 15:00
  • I needed to use lat_1 and lat_2 instead of lat1 and lat2. I suspect this applies post PROJ 6.0.0 – oortCloud Apr 26 at 1:58
18

It wasn't readily apparent to me how to use @sgillies answer, so here is a more verbose version:

import pyproj    
import shapely
import shapely.ops as ops
from shapely.geometry.polygon import Polygon
from functools import partial


geom = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
geom_area = ops.transform(
    partial(
        pyproj.transform,
        pyproj.Proj(init='EPSG:4326'),
        pyproj.Proj(
            proj='aea',
            lat1=geom.bounds[1],
            lat2=geom.bounds[3])),
    geom)

# Print the area in m^2
print geom_area.area 
  • 1
    The resulting value is not exactly the same as the one obtained in geojson.io. Why? – astrojuanlu Jul 6 '17 at 19:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.