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

Question

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})
Add coordinate reference system
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.

Answer 1

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

"""

# 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)



#------
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)