Geopandas – How to Determine a Polygon’s Area in Metric Units

areacoordinate systemgeopandasshape-area

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:

enter image description here

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:

enter image description here

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.

Best Answer

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)
Related Question