This question has come up several times, yet most solutions depend on a transformation to an equal area CRS. These solutions can be slow with large amounts of polygons, and also not accurate in high latitudes.
For example:
How to determine a polygon's area in a metric unit?
Area in KM from Polygon of coordinates
Getting polygon areas using GeoPandas
Calculating the Area by Square Feet with Geopandas
Is there a way to calculate the area of geographic coordinates without transforming?
Since version 0.7.0 geopandas has embedded the pyproj library as the crs object. pyproj, since version 2.3.0, has the ability to calculate the area of arbitrary polygons on a sphere. (see https://pyproj4.github.io/pyproj/stable/api/geod.html). The source of the math for this method is ultimately the geographiclib library. Thus, there is a straightforward way to do this calculation with minimal overhead and no new dependencies.
There is also an alternative version detailed in this answer which I use here for comparison. This is based on the line integral and Green's theorem. I've adapted it slightly to work with shapely polygons.
The following code block describes these two implementations. Both return the area, in meters^2, of polygons in geographic coordinates.
import numpy as np
import geopandas as gpd
from shapely.geometry.polygon import orient
def gpd_geographic_area(geodf):
if not geodf.crs and geodf.crs.is_geographic:
raise TypeError('geodataframe should have geographic coordinate system')
geod = geodf.crs.get_geod()
def area_calc(geom):
if geom.geom_type not in ['MultiPolygon','Polygon']:
return np.nan
# For MultiPolygon do each separately
if geom.geom_type=='MultiPolygon':
return np.sum([area_calc(p) for p in geom.geoms])
# orient to ensure a counter-clockwise traversal.
# See https://pyproj4.github.io/pyproj/stable/api/geod.html
# geometry_area_perimeter returns (area, perimeter)
return geod.geometry_area_perimeter(orient(geom, 1))[0]
return geodf.geometry.apply(area_calc)
from shapely.geometry import Polygon
def line_integral_polygon_area(geom, radius = 6378137):
"""
Computes area of spherical polygon, assuming spherical Earth.
Returns result in ratio of the sphere's area if the radius is specified.
Otherwise, in the units of provided radius.
lats and lons are in degrees.
from https://stackoverflow.com/a/61184491/6615512
"""
if geom.geom_type not in ['MultiPolygon','Polygon']:
return np.nan
# For MultiPolygon do each separately
if geom.geom_type=='MultiPolygon':
return np.sum([line_integral_polygon_area(p) for p in geom.geoms])
# parse out interior rings when present. These are "holes" in polygons.
if len(geom.interiors)>0:
interior_area = np.sum([line_integral_polygon_area(Polygon(g)) for g in geom.interiors])
geom = Polygon(geom.exterior)
else:
interior_area = 0
# Convert shapely polygon to a 2 column numpy array of lat/lon coordinates.
geom = np.array(geom.boundary.coords)
lats = np.deg2rad(geom[:,1])
lons = np.deg2rad(geom[:,0])
# Line integral based on Green's Theorem, assumes spherical Earth
#close polygon
if lats[0]!=lats[-1]:
lats = np.append(lats, lats[0])
lons = np.append(lons, lons[0])
#colatitudes relative to (0,0)
a = np.sin(lats/2)**2 + np.cos(lats)* np.sin(lons/2)**2
colat = 2*np.arctan2( np.sqrt(a), np.sqrt(1-a) )
#azimuths relative to (0,0)
az = np.arctan2(np.cos(lats) * np.sin(lons), np.sin(lats)) % (2*np.pi)
# Calculate diffs
# daz = np.diff(az) % (2*pi)
daz = np.diff(az)
daz = (daz + np.pi) % (2 * np.pi) - np.pi
deltas=np.diff(colat)/2
colat=colat[0:-1]+deltas
# Perform integral
integrands = (1-np.cos(colat)) * daz
# Integrate
area = abs(sum(integrands))/(4*np.pi)
area = min(area,1-area)
if radius is not None: #return in units of radius
return (area * 4*np.pi*radius**2) - interior_area
else: #return in ratio of sphere total area
return area - interior_area
# a wrapper to apply the method to a geo data.frame
def gpd_geographic_area_line_integral(geodf):
return geodf.geometry.apply(line_integral_polygon_area)
To compare the accuracy I used the Natural Earth States 10m shapefile (https://www.naturalearthdata.com/downloads/10m-cultural-vectors/). These more complex polygons are a better real world test as opposed to using squares or rectangles.
#------------
# import Natural Earth state polygons at 10m scale
#------------
import pooch
ne_states_link = 'https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_1_states_provinces.zip'
file_paths = pooch.retrieve(
ne_states_link,
known_hash=None,
processor=pooch.Unzip()
)
ne_states_shapefile = [p for p in file_paths if '.shp' in p][0]
All code blocks below assume the above two blocks are loaded.
I selected 15 different states/provinces around the world and obtained their "true" area from wikipedia. Note it is not clear how accurate this "true" area is, but it provides a decent baseline. I made sure to select some areas at high, low, and mid latitudes, and some with numerous islands.
#-------------------
# The "true" areas in km^2 of various states
#-------------------
state_total_land_areas = {
'Idaho':216443, # https://www.census.gov/geographies/reference-files/2010/geo/state-area.html
'Tennessee':109153,
'New Hampshire': 24214,
'Lapland':92674 , # https://en.wikipedia.org/wiki/Regions_of_Finland
'Kainuu':20197 ,
'Pirkanmaa': 12585,
'Maharashtra':307713, # https://en.wikipedia.org/wiki/States_and_union_territories_of_India
'Rajasthan': 342269,
'Goa': 3702,
'Otago': 31186, # https://en.wikipedia.org/wiki/Regions_of_New_Zealand
'Taranaki': 7254,
'Canterbury': 44504,
'Magallanes y Antártica Chilena': 132291.1, # https://en.wikipedia.org/wiki/Regions_of_Chile
'Aisén del General Carlos Ibáñez del Campo' : 108494.4,
'Los Lagos': 48583.6
}
state_total_land_areas = gpd.pd.DataFrame(state_total_land_areas.items(), columns=['name','true_area'])
# convert to m^2 for comparison
state_total_land_areas['true_area'] = state_total_land_areas['true_area'] * 1000**2
For the comparison I calculated percentage error for the 2 methods.
# Calculate the percent difference from "true" areas.
complex_polygons = gpd.read_file(ne_states_shapefile)
complex_polygons = complex_polygons[complex_polygons.name.isin(state_total_land_areas.name)]
complex_polygons = complex_polygons.merge(state_total_land_areas, on='name', how='left')
complex_polygons['area_pyproj'] = (gpd_geographic_area(complex_polygons) - complex_polygons.true_area) / complex_polygons.true_area
complex_polygons['area_line_integral'] = (gpd_geographic_area_line_integral(complex_polygons)- complex_polygons.true_area) / complex_polygons.true_area
print(complex_polygons[['name','area_pyproj','area_line_integral']])
Both methods get extremely close to truth, less than 1% in many cases. Neither has a lower error in all cases. The magnitude of error for each is approximately the same too. So one method is not more accurate than the other, they are essentially equal.
# The percent difference between calculated and "true" areas for the respective method.
name area_pyproj area_line_integral
0 Lapland 0.060416 0.055373
1 Kainuu 0.211758 0.206638
2 Idaho -0.000876 -0.000712
3 New Hampshire 0.002167 0.002488
4 Los Lagos -0.024568 -0.023877
5 Aisén -0.043047 -0.043341
6 Magallanes -0.079049 -0.080632
7 Rajasthan 0.002377 0.006424
8 Maharashtra -0.004749 0.000466
9 Goa -0.076008 -0.070663
10 Otago 0.055592 0.055534
11 Canterbury 0.030486 0.030845
12 Taranaki 0.104515 0.105967
13 Pirkanmaa 0.195929 0.191524
14 Tennessee 0.001336 0.003461
Despite both methods being the same, some errors are still larger than others. What causes the estimates to be off? I tested 3 potential causes.
def get_num_vertices(geom):
try:
if geom.geom_type=='MultiPolygon':
return np.sum([get_num_vertices(p) for p in geom.geoms])
return len(geom.boundary.coords)
except NotImplementedError:
# cannot get boundary when there is interior coordinates
return np.nan
complex_polygons['center_lat'] = [c.y for c in complex_polygons.centroid]
complex_polygons['num_vertices'] = complex_polygons.geometry.apply(get_num_vertices)
pivoted = complex_polygons.melt(
id_vars=['name','center_lat','num_vertices','true_area'],
value_vars=['area_pyproj','area_line_integral'],
var_name='method',
value_name='percent_error')
import seaborn as sns
sns.relplot(x="center_lat", y="percent_error",hue='method',data=pivoted).set(title='Percent error vs latitude')
sns.relplot(x="num_vertices", y="percent_error",hue='method', data=pivoted).set(title='Percent error vs num_vertices')
sns.relplot(x="true_area", y="percent_error",hue='method', data=pivoted).set(title='Percent error vs polygon area')
These results are with some figures.
Only latitude shows any kind of correlation with error. Neither larger sizes nor more points within the polygon contribute to higher error. Also note that the error is still extremely small here.
I did another test here where I recalculated areas for the full Natural Earth shapefile and compared there disagreement amongst all polygons.
complex_polygons = gpd.read_file(ne_states_shapefile)
complex_polygons['area_pyproj'] = gpd_geographic_area(complex_polygons) / 1000**2
complex_polygons['area_line_integral'] = gpd_geographic_area_line_integral(complex_polygons) / 1000**2
complex_polygons['method_difference'] = complex_polygons['area_pyproj'] - complex_polygons['area_line_integral']
complex_polygons['method_percent_difference'] = complex_polygons['method_difference'] / complex_polygons['area_pyproj']
complex_polygons['center_lat'] = [c.y for c in complex_polygons.centroid]
complex_polygons['num_vertices'] = complex_polygons.geometry.apply(get_num_vertices)
First off, Antarctica has the highest disagreement. Since this land mass is centered around a pole, its likely the calculus begins to fail there.
# Antarctica results
pyproj area: 12335385.988989946
line integral area:12263092.42375752
percent difference: 0.005860665024746868
Excluding Antarctica there is high agreement between the two methods.
Looking at percentage differences to highlight those small discrepancies some interesting patterns emerge. There is an interesting pattern of disagreement between them related to latitude. Around 45 latitude, north and south, is when they essentially agree. There is likely a mathematical explanation for that.
There is no correlation between the disagreement percent and area or number of vertices.
Here is the code to generate the figures.
# the rest minus antarctica
import seaborn as sns
complex_polygons = complex_polygons[complex_polygons.name!='Antarctica']
sns.relplot(x="area_pyproj", y="area_line_integral", data=complex_polygons).set(title='Both methods compared')
sns.relplot(x="center_lat", y="method_percent_difference", data=complex_polygons).set(title='Method differences vs latitude')
sns.relplot(x="num_vertices", y="method_percent_difference", data=complex_polygons).set(title='Method differences vs num_vertices')
sns.relplot(x="area_pyproj", y="method_percent_difference", data=complex_polygons).set(title='Method differences vs polgon area')
For a speed comparison I recalculated the full Natural Earth shapefile.
import timeit
complex_polygons = gpd.read_file(ne_states_shapefile)
pyproj_speed = timeit.timeit(lambda: gpd_geographic_area(complex_polygons), number=1)
line_integral_speed = timeit.timeit(lambda: gpd_geographic_area_line_integral(complex_polygons), number=1)
print('pyproj speed: {}'.format(round(pyproj_speed,2)))
print('line integral speed: {}'.format(round(line_integral_speed,2)))
The line integral method is clearly faster.
pyproj speed: 12.61
line integral speed: 1.53
One final thing to mention is instances where the polygons represent a regular grid of cells which are squares, and where the sides are latitudinal/longitudinal lines. I wrote about this in this answer, and will quickly compare it here for a speed test.
The method implemented
from math import radians, sin
def lat_lon_cell_area(lat_lon_grid_cell):
"""
Calculate the area of a cell, in meters^2, on a lat/lon grid.
This applies the following equation from Santinie et al. 2010.
S = (λ_2 - λ_1)(sinφ_2 - sinφ_1)R^2
S = surface area of cell on sphere
λ_1, λ_2, = bands of longitude in radians
φ_1, φ_2 = bands of latitude in radians
R = radius of the sphere
Santini, M., Taramelli, A., & Sorichetta, A. (2010). ASPHAA: A GIS‐Based
Algorithm to Calculate Cell Area on a Latitude‐Longitude (Geographic)
Regular Grid. Transactions in GIS, 14(3), 351-377.
https://doi.org/10.1111/j.1467-9671.2010.01200.x
Parameters
----------
lat_lon_grid_cell
A shapely box with coordinates on the lat/lon grid
Returns
-------
float
The cell area in meters^2
"""
# mean earth radius - https://en.wikipedia.org/wiki/Earth_radius#Mean_radius
AVG_EARTH_RADIUS_METERS = 6371008.8
west, south, east, north = lat_lon_grid_cell.bounds
west = radians(west)
east = radians(east)
south = radians(south)
north = radians(north)
return (east - west) * (sin(north) - sin(south)) * (AVG_EARTH_RADIUS_METERS**2)
Setting up a regular grid
from shapely.geometry import box
# Latitude (y)
lat_min, lat_max = -80, 80
# Longtidue (x)
lon_min, lon_max = -180, 180
# Edge length
side_length = 2
all_lats, all_lons = np.meshgrid(
np.arange(lat_min, lat_max, side_length),
np.arange(lon_min, lon_max, side_length)
)
polygons = []
for lon, lat in zip(all_lons.flatten(), all_lats.flatten()):
polygons.append(
box(lon, lat, lon+side_length, lat+side_length)
)
grid = gpd.GeoDataFrame({'geometry': gpd.GeoSeries(polygons)})
grid = grid.set_crs("EPSG:4326")
import timeit
cell_method = timeit.timeit(lambda: grid.geometry.apply(lat_lon_cell_area), number=10)
pyproj_speed = timeit.timeit(lambda: gpd_geographic_area(grid), number=10)
line_integral_speed = timeit.timeit(lambda: gpd_geographic_area_line_integral(grid), number=10)
print('pyproj regular grid speed: {}'.format(round(pyproj_speed,2)))
print('line integral regular speed: {}'.format(round(line_integral_speed,2)))
print('cell method regular grid speed: {}'.format(round(cell_method,2)))
The cell method, meant only for regular grids and not complex polygons, is significantly quicker than the other two methods. Interestingly, the magnitude of speed increase in the line integral over the pyproj method is not as high as it is with the more complex polygons. This suggests the speed up there is likely due to the vectorized methods used in the line integral method, which offer little speed improvement when polygons consist of only 4 points.
pyproj regular grid speed: 18.79
line integral regular speed: 16.35
cell method regular grid speed: 5.85
With the regular grid polygons I can also compare the agreement between methods.
grid['area_cell_method'] = grid.geometry.apply(lat_lon_cell_area) / 1000**2
grid['area_pyproj'] = gpd_geographic_area(grid) / 1000**2
grid['area_line_integral'] = gpd_geographic_area_line_integral(grid) / 1000**2
import seaborn as sns
sns.relplot(x="area_cell_method", y="area_pyproj", data=grid).set(title='cell vs pyproj method for regulard grid')
sns.relplot(x="area_cell_method", y="area_line_integral", data=grid).set(title='cell vs line integral method for regulard grid')
The cell method has near perfect agreement with the pyproj method. But, for the largest cells, has some large disagreement with the line integral method. Remember these are not country or province polygons but grid cells on a uniform lat/lon grid. The largest cells in a regular are at the equator, where the line integral method vastly overestimates some.
I traced this back to a grid cell with coords right at the equator. This is likely an edge condition for when grid cells have one of their bounds at 0.
# sort to get the highest value from the grid cell areas
grid.sort_values('area_line_integral').tail(2).geometry.values[0].bounds
(-2.0, -2.0, 0.0, 0.0)
If you want to calculate the area of polygons with geographic coordinates in geopandas:
lat_lon_cell_area) is extremely quick, straightforward, and produces approximately the same results as pyproj.As of Fall 2021:
Versions used here: