Tiling of the sphere by spherical triangles

Spherical Polygons and Areas¶

Overview¶

Determine the calculations of a spherical polygons based on a unit sphere.

Determine clockwise/counterclockwise ordering of points on spherical polygon
Area and Permieter of quadrilateral patch on a unit sphere
Determine if a given point is within a spherical polygon
Mean center of spherical polygon

Prerequisites¶

Concepts	Importance	Notes
Numpy	Necessary
Pandas	Necessary
Intro to Cartopy	Helpful	Will be used for plotting
Matplotlib	Helpful	Will be used for plotting

Time to learn: 40 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()

location_df.index = location_df["name"]

Determine clockwise/counterclockwise ordering of points on spherical polygon¶

True: when input points are in a clockwise order
False: when input points are in a counterclockwise (or co-linear) order

Shoelace Formula for Signed Polygon Area¶

TODO

def is_clockwise(pt_lst=None):
    # signed polygon area -> shoelace formula
    # positive = counterclockwise, negative = clockwise
    area = 0
    for i in range(0, len(pt_lst)):
        if i+1 < len(pt_lst):
            area += location_df.loc[pt_lst[i], "latitude"] * location_df.loc[pt_lst[i+1], "longitude"]
            area -= location_df.loc[pt_lst[i+1], "latitude"]  * location_df.loc[pt_lst[i], "longitude"]
        #area /= 2 # determine full sign area, unneeded when just working with signs
    if area < 0: 
        print("clockwise -> negative")
        return True
    if area > 0:
        print("counterclockwise -> positive")
        return False
    if area == 0:
        print("non-collinear -> zero") #https://en.wikipedia.org/wiki/Curve_orientation
        return False

is_clockwise(["boulder", "boston", "houston"])

clockwise -> negative

True

def plot_clockwise(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 arrow between points in order
    for i, pt in enumerate(pt_lst):
        if i+1 < len(pt_lst):
            ax.quiver(location_df.loc[pt_lst[i], "longitude"],
                      location_df.loc[pt_lst[i], "latitude"], 
                      (location_df.loc[pt_lst[i+1], "longitude"]-location_df.loc[pt_lst[i], "longitude"]), 
                      (location_df.loc[pt_lst[i+1], "latitude"]-location_df.loc[pt_lst[i], "latitude"]), 
                      angles='xy', scale_units='xy', scale=1)   
   # 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")
    if is_clockwise(pt_lst):
        clockwise = "Clockwise"
    else:
        clockwise = "Counterclockwise"
    plt.title(clockwise)
    plt.show()

plot_clockwise(["boulder", "boston", "houston"], -130, -60, 20, 60)

clockwise -> negative

/home/runner/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/IPython/core/formatters.py:402, in BaseFormatter.__call__(self, obj)
    400     pass
    401 else:
--> 402     return printer(obj)
    403 # Finally look for special method names
    404 method = get_real_method(obj, self.print_method)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/IPython/core/pylabtools.py:170, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    167     from matplotlib.backend_bases import FigureCanvasBase
    168     FigureCanvasBase(fig)
--> 170 fig.canvas.print_figure(bytes_io, **kw)
    171 data = bytes_io.getvalue()
    172 if fmt == 'svg':

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2154     # we do this instead of `self.figure.draw_without_rendering`
   2155     # so that we can inject the orientation
   2156     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157         self.figure.draw(renderer)
   2158 if bbox_inches:
   2159     if bbox_inches == "tight":

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     92 @wraps(draw)
     93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94     result = draw(artist, renderer, *args, **kwargs)
     95     if renderer._rasterizing:
     96         renderer.stop_rasterizing()

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/figure.py:3257, in Figure.draw(self, renderer)
   3254             # ValueError can occur when resizing a window.
   3256     self.patch.draw(renderer)
-> 3257     mimage._draw_list_compositing_images(
   3258         renderer, self, artists, self.suppressComposite)
   3260     renderer.close_group('figure')
   3261 finally:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    132 if not_composite or not has_images:
    133     for a in artists:
--> 134         a.draw(renderer)
    135 else:
    136     # Composite any adjacent images together
    137     image_group = []

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/mpl/geoaxes.py:509, in GeoAxes.draw(self, renderer, **kwargs)
    504         self.imshow(img, extent=extent, origin=origin,
    505                     transform=factory.crs, *factory_args[1:],
    506                     **factory_kwargs)
    507 self._done_img_factory = True
--> 509 return super().draw(renderer=renderer, **kwargs)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
   3223 if artists_rasterized:
   3224     _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
   3227     renderer, self, artists, self.get_figure(root=True).suppressComposite)
   3229 renderer.close_group('axes')
   3230 self.stale = False

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    132 if not_composite or not has_images:
    133     for a in artists:
--> 134         a.draw(renderer)
    135 else:
    136     # Composite any adjacent images together
    137     image_group = []

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/mpl/feature_artist.py:184, in FeatureArtist.draw(self, renderer)
    179     geoms = self._feature.geometries()
    180 else:
    181     # For efficiency on local maps with high resolution features (e.g
    182     # from Natural Earth), only create paths for geometries that are
    183     # in view.
--> 184     geoms = self._feature.intersecting_geometries(extent)
    186 stylised_paths = {}
    187 # Make an empty placeholder style dictionary for when styler is not
    188 # used.  Freeze it so that we can use it as a dict key.  We will need
    189 # to unfreeze all style dicts with dict(frozen) before passing to mpl.

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/feature/__init__.py:312, in NaturalEarthFeature.intersecting_geometries(self, extent)
    305 """
    306 Returns an iterator of shapely geometries that intersect with
    307 the given extent.
    308 The extent is assumed to be in the CRS of the feature.
    309 If extent is None, the method returns all geometries for this dataset.
    310 """
    311 self.scaler.scale_from_extent(extent)
--> 312 return super().intersecting_geometries(extent)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/feature/__init__.py:112, in Feature.intersecting_geometries(self, extent)
    109 if extent is not None and not np.isnan(extent[0]):
    110     extent_geom = sgeom.box(extent[0], extent[2],
    111                             extent[1], extent[3])
--> 112     return (geom for geom in self.geometries() if
    113             geom is not None and extent_geom.intersects(geom))
    114 else:
    115     return self.geometries()

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/feature/__init__.py:297, in NaturalEarthFeature.geometries(self)
    295     if reader.crs is not None:
    296         self._crs = reader.crs
--> 297     geometries = tuple(reader.geometries())
    298     _NATURAL_EARTH_GEOM_CACHE[key] = geometries
    299 else:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/shapereader.py:173, in BasicReader.geometries(self)
    161 def geometries(self):
    162     """
    163     Return an iterator of shapely geometries from the shapefile.
    164 
   (...)    171 
    172     """
--> 173     for shape in self._reader.iterShapes(bbox=self._bbox):
    174         # Skip the shape that can not be represented as geometry.
    175         if shape.shapeType != shapefile.NULL:
    176             yield sgeom.shape(shape)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:2806, in Reader.iterShapes(self, bbox)
   2802 if self.numShapes:
   2803     # Iterate exactly the number of shapes from shx header
   2804     for i in range(self.numShapes):
   2805         # MAYBE: check if more left of file or exit early?
-> 2806         shape = self.__shape(oid=i, bbox=bbox)
   2807         if shape:
   2808             yield shape

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:2684, in Reader.__shape(self, oid, bbox)
   2681 shapeType = unpack("<i", b_io.read(4))[0]
   2683 ShapeClass = SHAPE_CLASS_FROM_SHAPETYPE[shapeType]
-> 2684 shape = ShapeClass.from_byte_stream(
   2685     shapeType, b_io, recLength_bytes, oid=oid, bbox=bbox
   2686 )
   2688 # Seek to the end of this record as defined by the record header because
   2689 # the shapefile spec doesn't require the actual content to meet the header
   2690 # definition.  Probably allowed for lazy feature deletion.
   2691 # f.seek(next_shape)
   2693 return shape

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:1194, in _CanHaveBBox.from_byte_stream(cls, shapeType, b_io, next_shape, oid, bbox)
   1188 # else:
   1189 #     parts = None
   1190 #     partTypes = None
   1192 if nPoints:
   1193     kwargs["points"] = cast(
-> 1194         PointsT, cls._read_points_from_byte_stream(b_io, nPoints)
   1195     )
   1197     if shapeType in _HasZ_shapeTypes:
   1198         kwargs["zbox"], kwargs["z"] = _HasZ._read_zs_from_byte_stream(
   1199             b_io, nPoints
   1200         )

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:1136, in _CanHaveBBox._read_points_from_byte_stream(b_io, nPoints)
   1132 @staticmethod
   1133 def _read_points_from_byte_stream(
   1134     b_io: ReadableBinStream, nPoints: int
   1135 ) -> list[Point2D]:
-> 1136     flat = unpack(f"<{2 * nPoints}d", b_io.read(16 * nPoints))
   1137     return list(zip(*(iter(flat),) * 2))

error: unpack requires a buffer of 256 bytes

<Figure size 1500x1000 with 2 Axes>

plot_clockwise(["houston", "boston", "boulder"], -130, -60, 20, 60)

counterclockwise -> positive

plot_clockwise(["boulder", "boston", "greenwich", "cairo", "timbuktu"])

counterclockwise -> positive

Area and Perimeter of quadrilateral patch¶

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

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)

    area_ellps, perimeter_ellps = area_of_polygon_ellps(pt_lst)
    area_us, perimeter_us = area_of_polygon_unit_sphere(pt_lst)
    print(f"Ellipsoid Area   = {area_ellps} km^2")
    print(f"Unit 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()

area_ellps, perimeter_ellps = area_of_polygon_ellps(["boulder", "boston",
                                             "arecibo", "houston"])
area_us, perimeter_us = area_of_polygon_unit_sphere(["boulder", "boston",
                                             "arecibo", "houston"])
print(f"Area Ellipsoid   = {area_ellps} km^2")
print(f"Area Unit Sphere = {area_us} km^2")
print(f"Perimeter Ellipsoid = {perimeter_ellps} km")
print(f"Perimeter Unit 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 Unit Sphere = 5344606.94796931 km^2
Perimeter Ellipsoid = 10171.738963248145 km
Perimeter Unit SPhere = 10170.504728302833 km
Roughly 1.05% of the Earth's Surface
Roughly 1.05% of the Earth's Surface

plot_area(["boulder", "boston", "greenwich", "cairo", "arecibo", "houston"])

Ellipsoid Area   = 21872449.378265787 km^2
Unit Sphere Area = 21896220.663299154 km^2

plot_area(["redwoods", "rockford", "boston", "houston",], -130, -60, 20, 60)

Ellipsoid Area   = 3150946.426714995 km^2
Unit Sphere Area = 3149017.3086414044 km^2

plot_area(["redwoods", "boston", "houston"], -130, -60, 20, 60)

Ellipsoid Area   = 3788155.432965353 km^2
Unit Sphere Area = 3782548.632737316 km^2

TODO¶

Fix invalid overlapping polygon by re-ordering points into a clockwise order.

plot_area(["boulder", "boston", "houston", "boston", "cairo", "arecibo", "greenwich"])

Ellipsoid Area   = 914381.1786067598 km^2
Unit Sphere Area = 954445.989927043 km^2

Determine if a given point is within a spherical polygon¶

Single or list of points

def polygon_contains_points(pt_lst=None, polygon_pts=None, tolerance_m=1):
    # tolerance in meters
    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))
    polygon = Polygon(lat_lon_coords)
    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

def plot_polygon_pts(pt_lst=None, polygon_pts=None, tolerance_m=1,
                   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, tolerance_m)
    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 (tolerance {tolerance_m} m) = {pt_lst[contains_pts]}, not contained = {pt_lst[~contains_pts]}")
    plt.show()

polygon_contains_points(["boulder"], ["redwoods", "boston", "houston"], 1)

array([ True])

plot_polygon_pts(["boulder"], ["redwoods", "boston", "houston"], 1,
               -130, -60, 20, 60)

polygon_contains_points(["cape canaveral"], ["redwoods", "boston", "houston"], 1)

array([False])

plot_polygon_pts(["cape canaveral"], ["redwoods", "boston", "houston"], 1,
               -130, -60, 20, 60)

plot_polygon_pts(["boulder", "cape canaveral"], ["redwoods", "boston", "houston"], 1,
               -130, -60, 20, 60)

plot_polygon_pts(["boulder", "redwoods"], ["rockford", "boston", "cape canaveral"], 1,
               -130, -60, 20, 60)

Mean center of spherical polygon¶

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)

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)

plot_centroid(["redwoods", "boulder", "cape canaveral", "houston"],
               -130, -60, 20, 60)

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