How to determine a polygon's area in a metric unit?

I have created a square grid and I would like to figure out the real world area of each polygon in a metric unit e.q. square metres via `GeoDataFrame.area`.

``````import geopandas as gpd
from shapely.geometry import Polygon

# Latitude (y)
lat_min = -90
lat_max = 90

# Longtidue (x)
lon_min = -180
lon_max = 180

# Edge length
side_length = 0.5

# List to which we will append the polygons
temp_polygons = []

for y in range(int(abs(lat_min-lat_max)/side_length)):

for x in range(int(abs(lon_min-lon_max)/side_length)):

x_start = lon_min + x * side_length
y_start = lat_min + y * side_length

bottom_left_x = x_start
bottom_left_y = y_start

bottom_right_x = x_start + side_length
bottom_right_y = y_start

top_right_x = x_start + side_length
top_right_y = y_start + side_length

top_left_x = x_start
top_left_y = y_start + side_length

# Append polygon to list
temp_polygons.append(
# Polygon
Polygon([
# Bottom-left
(bottom_left_x,bottom_left_y),
# Bottom-right
(bottom_right_x,bottom_right_y),
# Top-left
(top_right_x,top_right_y),
# Top-right
(top_left_x,top_left_y)
])
)

# Turn list of polygons into GeoSeries
temp_polygons = gpd.GeoSeries(temp_polygons)

# Turn GeoSeries into GeoDataFrame
raster = gpd.GeoDataFrame({'geometry': temp_polygons})

raster.crs = {"init":"epsg:4326"}

# Change projection
# https://en.wikipedia.org/wiki/Web_Mercator_projection#EPSG:3857
raster['geometry'] = (
raster['geometry'].to_crs({'init': 'EPSG:3857'})
)

# Get area
raster['area'] = raster.area

``````

First, the geometries are as follows:

However, after applying ...

``````# Change projection
# https://en.wikipedia.org/wiki/Web_Mercator_projection#EPSG:3857
raster['geometry'] = (
raster['geometry'].to_crs({'init': 'EPSG:3857'})
)
``````

... there are a lot of inf values and hence the cells' area is `nan`:

What am I doing wrong? Did I select a wrong CRS? How can I determine the area of each grid cell in a metric unit?

Edit:

This answer seems to work for me. However, with ~25 iterations/second it's rather slow (~3 hours for ~260,000 rows). It's not prohibitively slow but there must be a quicker solution.

• Thank you for your comment. This answer seems to work: gis.stackexchange.com/a/166421/97137; However, it's rather slow. ~25 iterations per second. That's too slow for my >250,000 rows. Oct 5, 2021 at 11:41
• Tangential comment: Do you remember why you thought that EPSG:3857 would give you actual meters? I would love to fix the web on that misleading thought so if it was a website or something, I'd be interested. Oct 5, 2021 at 16:34
• While EPSG:3857 won't give you useful areal measurements, are you sure your axes are correct? I get no NaNs with pyproj==3.1.0 and geopandas==0.9.0. Oct 5, 2021 at 20:03
• @bugmenot123 Thank you for your comments! Sorry but I cannot recall where I stumbled upon that information. It may have even been Stackexchange but I cannot find it anymore. Regarding the 2nd comment I am not a hundred percent sure what you mean. I provided all the code in the minimal working example. Oct 5, 2021 at 20:15

Here is a way to calculate lat/long cell area in m^2 without any transformation as long as they are regular grid cells. In other words the polygons must be squares where each side is a latitudinal or longitudinal line. From the article Santini et al. 2010

It's quite fast. It did the ~260,000 cells in your `raster` data in ~10 seconds on my small desktop computer.

``````timeit.timeit(lambda: raster.geometry.apply(lat_lon_cell_area), number=1)
10.299032560084015

raster.shape
(259200, 2)
``````

I confirmed the values by writing the `raster` grid to a shapefile and then opening it in qgis where I measured some cell areas manually. The two were approximately the same.

``````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 Santini 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

"""

west, south, east, north = lat_lon_grid_cell.bounds

return (east - west) * (sin(north) - sin(south)) * (AVG_EARTH_RADIUS_METERS**2)

#------
# Example using a 2 degree lat/lon grid
#------

import numpy as np
import geopandas as gpd
from shapely.geometry import box

# Latitude (y)
lat_min = -90
lat_max = 90

# Longitude (x)
lon_min = -180
lon_max = 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)
)

raster = gpd.GeoDataFrame({'geometry': gpd.GeoSeries(polygons)})

raster['cell_area'] = raster.geometry.apply(lat_lon_cell_area)

``````
• Thank you for your superb answer! It seems to work perfectly fine. I checked `raster['cell_area'].sum()` and it deviates by only 1% (`510063456035226.2 / 510065880972871.75 = 0.9999952458344382`) from the result I cross-linked in my answer (I guess the deviation is caused by e.g. post-digit differences in e.g. the earth radius). But your approach is much much quicker! ... perfect answer! Thanks! Oct 6, 2021 at 6:26