How to calculate the area of a specific region defined by a set of points with latitudes and longitudes using Python ?

Introduction

In this document, we will discuss how to calculate the area of a specific region defined by a set of points with latitudes and longitudes using Python. This is a commonly used technique in various fields such as geography, geology, and remote sensing.

Define a geographical area by specifying its longitude and latitude coordinates

Before we dive into calculating the area, let's first understand what latitudes and longitudes are and how they help define a specific location on the Earth's surface. Latitudes and longitudes are imaginary lines that crisscross the globe, creating a grid-like structure. Latitudes run horizontally around the Earth while longitudes run vertically.

The combination of latitudes and longitudes allows us to pinpoint any location on the Earth's surface. This is extremely useful when calculating the area of a specific region as it helps define the boundaries of that region accurately.

Using Python, let's begin by defining a geographical area through the specification of its longitude and latitude coordinates:

    coords  = ((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0) ) # point: (longitude,latitude)

Calculating Area using Python

There are various methods to calculate the area of a specific region with python

Using shapely and pyproj

from pyproj import Geod
from shapely.geometry import Polygon

geod = Geod(ellps="WGS84")

coords  = ((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0) ) # point: (longitude,latitude)

poly = Polygon(coords)

area = abs(geod.geometry_area_perimeter(poly)[0])

print('Area {:.3f} m^2'.format(area))
print('Area {:.3f} km^2'.format(area/1000**2.0))

which returns

Area 12308778361.469 m^2
Area 12308.778 km^2

Since our area is approximately a square, we can calculate the length of a side using the following formula

import numpy as np

print( np.sqrt(area) / 1000 )

gives

110.94493391529625

in km.This makes sense because approximately 1 degree corresponds to 111.1 km at the equator (see How big is a degree ?).

Using shapely transform function

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

geom = Polygon([(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)])

geom_area = ops.transform(
    partial(
        pyproj.transform,
        pyproj.Proj(init='EPSG:4326'),
        pyproj.Proj(
            proj='aea',
            lat_1=geom.bounds[1],
            lat_2=geom.bounds[3]
        )
    ),
    geom)


area_m2 = geom_area.area

# Print the area in m^2
print(area_m2)

print( np.sqrt(area_m2) / 1000 )

returns

12308463846.396065
110.9435164684988

Using wkt loads to get area of a specific region

from pyproj import Geod
from shapely import wkt

# specify a named ellipsoid
geod = Geod(ellps="WGS84")

poly = wkt.loads('''\
POLYGON ((0.0 0.0, 0.0 1.0, 1.0 1.0, 1.0 0.0, 0.0 0.0))''')

area = abs(geod.geometry_area_perimeter(poly)[0])

print('Area: {:.3f} m^2'.format(area))

returns

Area: 12308778361.469 m^2

References