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

Links | Site |
---|---|

shapely.Polygon | shapely.readthedocs.io |

geopandas.GeoSeries.area | geopandas.org |

How big is a degree ? | sco.wisc.edu |