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:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

---------------------------------------------------------------------------
ConnectionResetError                      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:2155, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2152     # we do this instead of `self.figure.draw_without_rendering`
   2153     # so that we can inject the orientation
   2154     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2155         self.figure.draw(renderer)
   2156 if bbox_inches:
   2157     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:524, in GeoAxes.draw(self, renderer, **kwargs)
    519         self.imshow(img, extent=extent, origin=origin,
    520                     transform=factory.crs, *factory_args[1:],
    521                     **factory_kwargs)
    522 self._done_img_factory = True
--> 524 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:3216, in _AxesBase.draw(self, renderer)
   3213 if artists_rasterized:
   3214     _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3216 mimage._draw_list_compositing_images(
   3217     renderer, self, artists, self.get_figure(root=True).suppressComposite)
   3219 renderer.close_group('axes')
   3220 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:185, in FeatureArtist.draw(self, renderer)
    180     geoms = self._feature.geometries()
    181 else:
    182     # For efficiency on local maps with high resolution features (e.g
    183     # from Natural Earth), only create paths for geometries that are
    184     # in view.
--> 185     geoms = self._feature.intersecting_geometries(extent)
    187 stylised_paths = {}
    188 # Make an empty placeholder style dictionary for when styler is not
    189 # used.  Freeze it so that we can use it as a dict key.  We will need
    190 # 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:309, in NaturalEarthFeature.intersecting_geometries(self, extent)
    302 """
    303 Returns an iterator of shapely geometries that intersect with
    304 the given extent.
    305 The extent is assumed to be in the CRS of the feature.
    306 If extent is None, the method returns all geometries for this dataset.
    307 """
    308 self.scaler.scale_from_extent(extent)
--> 309 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:291, in NaturalEarthFeature.geometries(self)
    289 key = (self.name, self.category, self.scale)
    290 if key not in _NATURAL_EARTH_GEOM_CACHE:
--> 291     path = shapereader.natural_earth(resolution=self.scale,
    292                                      category=self.category,
    293                                      name=self.name)
    294     geometries = tuple(shapereader.Reader(path).geometries())
    295     _NATURAL_EARTH_GEOM_CACHE[key] = geometries

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/shapereader.py:306, in natural_earth(resolution, category, name)
    302 ne_downloader = Downloader.from_config(('shapefiles', 'natural_earth',
    303                                         resolution, category, name))
    304 format_dict = {'config': config, 'category': category,
    305                'name': name, 'resolution': resolution}
--> 306 return ne_downloader.path(format_dict)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/__init__.py:203, in Downloader.path(self, format_dict)
    200     result_path = target_path
    201 else:
    202     # we need to download the file
--> 203     result_path = self.acquire_resource(target_path, format_dict)
    205 return result_path

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/shapereader.py:361, in NEShpDownloader.acquire_resource(self, target_path, format_dict)
    357 url = self.url(format_dict)
    359 shapefile_online = self._urlopen(url)
--> 361 zfh = ZipFile(io.BytesIO(shapefile_online.read()), 'r')
    363 for member_path in self.zip_file_contents(format_dict):
    364     member = zfh.getinfo(member_path.replace('\\', '/'))

File ~/micromamba/envs/cookbook-gc/lib/python3.13/http/client.py:495, in HTTPResponse.read(self, amt)
    493 else:
    494     try:
--> 495         s = self._safe_read(self.length)
    496     except IncompleteRead:
    497         self._close_conn()

File ~/micromamba/envs/cookbook-gc/lib/python3.13/http/client.py:642, in HTTPResponse._safe_read(self, amt)
    635 def _safe_read(self, amt):
    636     """Read the number of bytes requested.
    637 
    638     This function should be used when <amt> bytes "should" be present for
    639     reading. If the bytes are truly not available (due to EOF), then the
    640     IncompleteRead exception can be used to detect the problem.
    641     """
--> 642     data = self.fp.read(amt)
    643     if len(data) < amt:
    644         raise IncompleteRead(data, amt-len(data))

File ~/micromamba/envs/cookbook-gc/lib/python3.13/socket.py:719, in SocketIO.readinto(self, b)
    717     raise OSError("cannot read from timed out object")
    718 try:
--> 719     return self._sock.recv_into(b)
    720 except timeout:
    721     self._timeout_occurred = True

File ~/micromamba/envs/cookbook-gc/lib/python3.13/ssl.py:1304, in SSLSocket.recv_into(self, buffer, nbytes, flags)
   1300     if flags != 0:
   1301         raise ValueError(
   1302           "non-zero flags not allowed in calls to recv_into() on %s" %
   1303           self.__class__)
-> 1304     return self.read(nbytes, buffer)
   1305 else:
   1306     return super().recv_into(buffer, nbytes, flags)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/ssl.py:1138, in SSLSocket.read(self, len, buffer)
   1136 try:
   1137     if buffer is not None:
-> 1138         return self._sslobj.read(len, buffer)
   1139     else:
   1140         return self._sslobj.read(len)

ConnectionResetError: [Errno 104] Connection reset by peer

<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