Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Spherical Polygons and Areas

Authors
Affiliations
NSF National Center for Atmospheric Research
Tiling of the sphere by spherical triangles

Spherical Polygons and Areas

Overview

In this notebook, we will calculate how to use spherical polygons on a sphere to determine areas and perimeters.

  • Calculate area and permieter of quadrilateral patch on a sphere

  • Determine if a given point is within a spherical polygon

  • Mean center of spherical polygon

Prerequisites

ConceptsImportanceNotes
NumpyNecessaryUsed to work with large arrays
PandasNecessaryUsed to read in and organize data (in particular dataframes)
Intro to CartopyHelpfulWill be used for adding maps to plotting
MatplotlibHelpfulWill be used for plotting

  • Time to learn: 20 minutes

Imports

  • Import Packages

  • Setup location dataframe with coordinates

import pandas as pd       # reading in data for location information from text file
import numpy as np        # working with arrays, vectors, cross/dot products, and radians

from pyproj import Geod   # working with the Earth as an ellipsod (WGS-84)

from shapely.geometry import Point
from shapely.geometry.polygon import Polygon

import matplotlib.pyplot as plt                        # plotting a graph
from cartopy import crs as ccrs, feature as cfeature   # plotting a world map
# Get all Coordinates for Locations
location_df = pd.read_csv("../location_full_coords.txt")
location_df = location_df.rename(columns=lambda x: x.strip()) # strip excess white space from column names and values
location_df.head()
Loading...
location_df.index = location_df["name"]

Area and Perimeter of quadrilateral patch

We can make use of the pyproj Python package to calculate the area and perimeter of a patch formed by a list of latitude and longitude coordinates for an ellipsoid or a sphere.

def area_of_polygon_ellps(poly_pts=None):
    geod = Geod(ellps="WGS84")
    longitudes = [location_df.loc[pt, "longitude"] for pt in poly_pts]
    latitudes = [location_df.loc[pt, "latitude"] for pt in poly_pts]
    poly_area_m, poly_perimeter_m = geod.polygon_area_perimeter(longitudes, latitudes)
    return abs(poly_area_m) * 1e-6, poly_perimeter_m/1000 # km^2 and km

def area_of_polygon_sphere(poly_pts=None):
    geod = Geod(ellps="sphere") # 'sphere': {'a': 6370997.0, 'b': 6370997.0, 'description': 'Normal Sphere (r=6370997)'
    longitudes = [location_df.loc[pt, "longitude"] for pt in poly_pts]
    latitudes = [location_df.loc[pt, "latitude"] for pt in poly_pts]
    poly_area_m, poly_perimeter_m = geod.polygon_area_perimeter(longitudes, latitudes)
    return abs(poly_area_m) * 1e-6, poly_perimeter_m/1000 # km^2 and km

Let’s compare how area and perimeter vary based on whether you use an ellpsoid like WGS-84 and a sphere. The changes will become more pronounced as you move further towards the poles.

area_ellps, perimeter_ellps = area_of_polygon_ellps(["boulder", "boston",
                                             "arecibo", "houston"])
area_us, perimeter_us = area_of_polygon_sphere(["boulder", "boston",
                                             "arecibo", "houston"])
print(f"Area Ellipsoid   = {area_ellps} km^2")
print(f"Area Sphere = {area_us} km^2")
print(f"Perimeter Ellipsoid = {perimeter_ellps} km")
print(f"Perimeter Sphere = {perimeter_us} km")
print(f"Roughly {(area_ellps/509600000)*100:.2f}% of the Earth's Surface")
print(f"Roughly {(area_us/509600000)*100:.2f}% of the Earth's Surface")
Area Ellipsoid   = 5342585.6476998255 km^2
Area Sphere = 5344606.94796931 km^2
Perimeter Ellipsoid = 10171.738963248145 km
Perimeter Sphere = 10170.504728302833 km
Roughly 1.05% of the Earth's Surface
Roughly 1.05% of the Earth's Surface

Plot Area of Quadrilateral Patch

Let’s see what spherical polygons looks like on a world map!

def plot_area(pt_lst=None,
                   lon_west=-180, lon_east=180,
                   lat_south=-90, lat_north=90):
    # Set up world map plot
    fig = plt.subplots(figsize=(15, 10))
    projection_map = ccrs.PlateCarree()
    ax = plt.axes(projection=projection_map)
    ax.set_extent([lon_west, lon_east, lat_south, lat_north], crs=projection_map)
    ax.coastlines(color="black")
    ax.add_feature(cfeature.STATES, edgecolor="black")

   # plot points
    longitudes = [location_df.loc[x, "longitude"] for x in pt_lst] # longitude
    latitudes = [location_df.loc[y, "latitude"] for y in pt_lst] # latitude
    plt.scatter(longitudes, latitudes, s=100, c="red")
    plt.fill(longitudes, latitudes, alpha=0.5)

    # determine the area and perimeter
    area_ellps, perimeter_ellps = area_of_polygon_ellps(pt_lst)
    area_us, perimeter_us = area_of_polygon_sphere(pt_lst)
    print(f"Ellipsoid Area   = {area_ellps} km^2")
    print(f"Sphere Area = {area_us} km^2")
    plt.title(f"Roughly {(area_ellps/509600000)*100:.2f}% ({(area_us/509600000)*100:.2f}%) of the Earth's Surface")
    plt.show()
plot_area(["boulder", "boston", "greenwich", "cairo", "arecibo", "houston"])
Ellipsoid Area   = 21872449.378265787 km^2
Sphere Area = 21896220.663299154 km^2
<Figure size 1500x1000 with 2 Axes>
plot_area(["redwoods", "rockford", "boston", "houston",], -130, -60, 20, 60)
Ellipsoid Area   = 3150946.426714995 km^2
Sphere Area = 3149017.3086414044 km^2
<Figure size 1500x1000 with 2 Axes>
plot_area(["redwoods", "boston", "houston"], -130, -60, 20, 60)
Ellipsoid Area   = 3788155.432965353 km^2
Sphere Area = 3782548.632737316 km^2
<Figure size 1500x1000 with 2 Axes>
# For example: here is a list of unordered points
plot_area(["boulder", "boston", "houston", "boston", "cairo", "arecibo", "greenwich"])
Ellipsoid Area   = 914381.1786067598 km^2
Sphere Area = 954445.989927043 km^2
<Figure size 1500x1000 with 2 Axes>

Determine if a given point is within a spherical polygon

Now that we have a spherical polygon, how can be determine if a new point lies within the area of the polygon? The shapely Python package makes this quite easy, let’s give it a shot!

def polygon_contains_points(pt_lst=None, polygon_pts=None):
    longitudes = [location_df.loc[pt, "longitude"] for pt in polygon_pts]
    latitudes = [location_df.loc[pt, "latitude"] for pt in polygon_pts]
    lat_lon_coords = tuple(zip(longitudes, latitudes))
    # setup polygon
    polygon = Polygon(lat_lon_coords)
    # check if new point lies within the polygon
    contains = np.vectorize(lambda pt: polygon.contains(Point((location_df.loc[pt, "longitude"],
                                                               location_df.loc[pt, "latitude"]))))
    contained_by_polygon = contains(np.array(pt_lst))
    return contained_by_polygon
polygon_contains_points(["boulder"], ["redwoods", "boston", "houston"])
array([ True])

Plot and See if New Point within Polygon

On a map, it can be fairly intuitive to see if a point lies within a polygon or not. Let’s give it a shot!

def plot_polygon_pts(pt_lst=None, polygon_pts=None,
                   lon_west=-180, lon_east=180,
                   lat_south=-90, lat_north=90):
    # Set up world map plot
    fig = plt.subplots(figsize=(15, 10))
    projection_map = ccrs.PlateCarree()
    ax = plt.axes(projection=projection_map)
    ax.set_extent([lon_west, lon_east, lat_south, lat_north], crs=projection_map)
    ax.coastlines(color="black")
    ax.add_feature(cfeature.STATES, edgecolor="black")

    # plot polygon points
    longitudes = [location_df.loc[x, "longitude"] for x in polygon_pts] # longitude
    latitudes = [location_df.loc[y, "latitude"] for y in polygon_pts] # latitude
    plt.scatter(longitudes, latitudes, s=50, c="blue")
    plt.fill(longitudes, latitudes, alpha=0.5)

    # plot check points
    pt_lst = np.array(pt_lst)
    contains_pts = polygon_contains_points(pt_lst, polygon_pts)
    longitudes = [location_df.loc[x, "longitude"] for x in pt_lst[contains_pts]] # longitude
    latitudes = [location_df.loc[y, "latitude"] for y in pt_lst[contains_pts]] # latitude
    plt.scatter(longitudes, latitudes, s=100, c="green", label="Within Polygon")
    longitudes = [location_df.loc[x, "longitude"] for x in pt_lst[~contains_pts]] # longitude
    latitudes = [location_df.loc[y, "latitude"] for y in pt_lst[~contains_pts]] # latitude
    plt.scatter(longitudes, latitudes, s=100, c="red", label="Not within Polygon")

    plt.legend(loc="lower left")
    plt.title(f"Points contained within polygon = {pt_lst[contains_pts]}, not contained = {pt_lst[~contains_pts]}")
    plt.show()
polygon_contains_points(["boulder"], ["redwoods", "boston", "houston"])
array([ True])
plot_polygon_pts(["boulder"], ["redwoods", "boston", "houston"],
               -130, -60, 20, 60)
<Figure size 1500x1000 with 2 Axes>
polygon_contains_points(["cape canaveral"], ["redwoods", "boston", "houston"])
array([False])
plot_polygon_pts(["cape canaveral"], ["redwoods", "boston", "houston"],
               -130, -60, 20, 60)
<Figure size 1500x1000 with 2 Axes>
plot_polygon_pts(["boulder", "cape canaveral"], ["redwoods", "boston", "houston"],
               -130, -60, 20, 60)
<Figure size 1500x1000 with 2 Axes>
plot_polygon_pts(["boulder", "redwoods"], ["rockford", "boston", "cape canaveral"],
               -130, -60, 20, 60)
<Figure size 1500x1000 with 2 Axes>

Mean center of spherical polygon

A spherical polygon can have a fairly regular shape, especially as more points are added. But it is fairly simple using the shapley Python package to determine the mean center.

def polygon_centroid(polygon_pts=None):
    longitudes = [location_df.loc[x, "longitude"] for x in polygon_pts]
    latitudes = [location_df.loc[y, "latitude"] for y in polygon_pts]
    lat_lon_coords = tuple(zip(longitudes, latitudes))
    polygon = Polygon(lat_lon_coords)
    return (polygon.centroid.y, polygon.centroid.x)
polygon_centroid(["boulder", "boston", "houston"])
(37.30896666666666, -90.47586666666665)

Plot Centroid

The center of a polygon can be fairly apparently on a map, let’s give it a look!

def plot_centroid(polygon_pts=None,
                   lon_west=-180, lon_east=180,
                   lat_south=-90, lat_north=90):
    # Set up world map plot
    fig = plt.subplots(figsize=(15, 10))
    projection_map = ccrs.PlateCarree()
    ax = plt.axes(projection=projection_map)
    ax.set_extent([lon_west, lon_east, lat_south, lat_north], crs=projection_map)
    ax.coastlines(color="black")
    ax.add_feature(cfeature.STATES, edgecolor="black")

   # plot polygon points
    longitudes = [location_df.loc[x, "longitude"] for x in polygon_pts] # longitude
    latitudes = [location_df.loc[y, "latitude"] for y in polygon_pts] # latitude
    plt.scatter(longitudes, latitudes, s=50, c="blue")
    plt.fill(longitudes, latitudes, alpha=0.5)

    # plot check point
    centeroid = polygon_centroid(polygon_pts)
    plt.scatter(centeroid[1], centeroid[0], s=100, c="red")
    plt.title(f"Centroid = {centeroid}")
    plt.show()
plot_centroid(["boulder", "boston", "houston"],
               -130, -60, 20, 60)
<Figure size 1500x1000 with 2 Axes>
plot_centroid(["redwoods", "boulder", "cape canaveral", "houston"],
               -130, -60, 20, 60)
<Figure size 1500x1000 with 2 Axes>

Summary

This notebook covers working with spherical polygons to determine the ordering of coordinates, center of polygons, and whether or not a point lies within a spherical polygon

Resources and references

Tutorials
Angles and Great Circles