Polygon Simplification#

This guide demonstrates customizing polygon requests to simplify the geometry of the returned polygons.

Contrails API supports both GeoJSON and KML polygon formats. We only work with GeoJSON polygons in this notebook.

This notebook uses Shapely to interact with polygons. In particular, GEOS is required to convert GeoJSON to Shapely polygons.

[1]:
import json
import os
from pprint import pprint

import matplotlib.lines as mlines
import matplotlib.pyplot as plt
import pandas as pd
import requests
import shapely
[2]:
# Define credentials
URL = os.environ["API_URL_BASE"]  # https://api.contrails.org/v0
api_key = os.environ["CONTRAILS_API_KEY"]
headers = {"x-api-key": api_key}

Query parameters#

The following query parameters are used in the /grid/cocip endpoint:

  • time: The time for the polygon prediction. Must be an ISO 8601 datetime string or unix timestamp (seconds since epoch).

  • format: Set format=geojson to return polygons in GeoJSON format.

  • aircraft_type: Specify an aircraft type for the polygon prediction.

  • threshold: The energy forcing threshold for the polygon boundary. See our reference for more information.

  • interiors: Set interiors=true to return interior polygons.

  • simplify: Used to control the extent to which polygon geometry is simplified.

  • convex_hull: Set convex_hull=true to return the convex hull of the polygon.

  • bbox: Used to downselect horizontally.

  • flight_level: Used to downselect vertically.

Default parameters#

Here, we request CoCiP polygons at a single flight level using the default polygon simplify parameters.

The flight_level query parameter was added in version 0.14.0.

[3]:
params = {
    "time": "2023-02-10",
    "aircraft_type": "A320",
    "threshold": 1e8,
    "format": "geojson",
    "flight_level": 360,
}

r = requests.get(f"{URL}/grid/cocip", params=params, headers=headers)
print(f"HTTP Response Code: {r.status_code} {r.reason}")
HTTP Response Code: 200 OK
[4]:
geojson = r.json()
feature = geojson["features"][0]

# Print out the metadata
pprint(feature["properties"])

# Convert the output to shapely
# Unfortunately, GEOS can't handle the altitude coordinate returned by the API
# Manually loop through and remove it
for coords in feature["geometry"]["coordinates"]:
    for poly in coords:
        for point in poly:
            del point[2]

multipoly = shapely.from_geojson(json.dumps(feature))
{'aircraft_type': 'A320',
 'cocip_dt_integration': '5 minutes',
 'cocip_max_contrail_age': '12 hours',
 'humidity_scaling_formula': 'rhi -> (rhi / rhi_adj) ^ rhi_boost_exponent',
 'humidity_scaling_name': 'exponential_boost_latitude_customization',
 'level': 360,
 'level_long_name': 'Flight Level',
 'level_standard_name': 'FL',
 'level_units': 'hectofeet',
 'long_name': 'Energy forcing per meter of flight trajectory',
 'met_source_dataset': 'ERA5',
 'met_source_product': 'reanalysis',
 'met_source_provider': 'ECMWF',
 'name': 'ef_per_m',
 'polygon_iso_value': 100000000.0,
 'pycontrails_version': '0.32.2',
 'time': '2023-02-10T00:00:00Z',
 'units': 'J / m'}

Polygon topology#

GeoJSON polygons have both interior and exterior boundaries. Below use red for exterior boundaries and blue for interior boundaries.

[5]:
for poly in multipoly.geoms:
    plt.plot(*poly.exterior.xy, color="red", lw=1)
    for interior in poly.interiors:
        plt.plot(*interior.xy, color="blue", lw=1)
../_images/notebooks_polygons_8_0.png
[6]:
# Print some statistics on the polygons
n_exterior_rings = len(multipoly.geoms)
n_interior_rings = sum(len(p.interiors) for p in multipoly.geoms)
n_exterior_vertices = sum(len(p.exterior.coords) for p in multipoly.geoms)
n_interior_vertices = sum(len(i.coords) for p in multipoly.geoms for i in p.interiors)
area = multipoly.area

print(f"Number of exterior rings: {n_exterior_rings}")
print(f"Number of interior rings: {n_interior_rings}")
print(f"Number of exterior vertices: {n_exterior_vertices}")
print(f"Number of interior vertices: {n_interior_vertices}")
print(f"Total area: {area}")  # units are in the image of the Plate carree projection
Number of exterior rings: 127
Number of interior rings: 10
Number of exterior vertices: 4389
Number of interior vertices: 117
Total area: 1481.7415000000037
[7]:
# Formalize some of cells above into functions for later use


def get_shapely(**kwargs):
    params = {
        "time": "2023-02-10",
        "aircraft_type": "A320",
        "threshold": 1e8,
        "format": "geojson",
        "flight_level": 360,
    }
    params.update(kwargs)

    r = requests.get(f"{URL}/grid/cocip", params=params, headers=headers)
    r.raise_for_status()

    geojson = r.json()
    (feature,) = geojson["features"]

    for coords in feature["geometry"]["coordinates"]:
        for poly in coords:
            for point in poly:
                del point[2]

    return shapely.from_geojson(json.dumps(feature))


def polygon_metrics(multipoly):
    n_exterior_rings = len(multipoly.geoms)
    n_interior_rings = sum(len(p.interiors) for p in multipoly.geoms)
    n_exterior_vertices = sum(len(p.exterior.coords) for p in multipoly.geoms)
    n_interior_vertices = sum(len(i.coords) for p in multipoly.geoms for i in p.interiors)
    area = multipoly.area
    return {
        "n_exterior_rings": n_exterior_rings,
        "n_interior_rings": n_interior_rings,
        "n_exterior_vertices": n_exterior_vertices,
        "n_interior_vertices": n_interior_vertices,
        "area": area,
    }

Interior polygons#

Set the query parameter interiors=false to only return exterior polygons.

As seen in the polygon_metrics output and in the plot, the exterior polygons remain the same as before, but the interior polygons are removed.

[8]:
multipoly = get_shapely(interiors=False)

print(polygon_metrics(multipoly))

for poly in multipoly.geoms:
    plt.plot(*poly.exterior.xy, color="red", lw=1)
    assert len(poly.interiors) == 0
{'n_exterior_rings': 127, 'n_interior_rings': 0, 'n_exterior_vertices': 4389, 'n_interior_vertices': 0, 'area': 1498.7391000000036}
../_images/notebooks_polygons_12_1.png

The simplify query parameter#

Use the simplify query parameter to customize the extent to which polygon boundaries are simplified. This parameter is an integer between 0 and 10, inclusive. A value of 0 means no simplification. A value of 10 simplifies aggressively (thereby introducing some error). The default value is 3. A higher value also drops more small polygons.

[9]:
multipoly_list = [get_shapely(interiors=False, simplify=s) for s in range(0, 11)]
[10]:
# Plot the least and most aggressively simplified polygons side by side

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))

ax1.set_title("simplify=0")
multipoly = multipoly_list[0]
for poly in multipoly.geoms:
    ax1.plot(*poly.exterior.xy, color="red", lw=1)

ax2.set_title("simplify=10")
multipoly = multipoly_list[-1]
for poly in multipoly.geoms:
    ax2.plot(*poly.exterior.xy, color="red", lw=1)
../_images/notebooks_polygons_15_0.png

Metrics#

Make a table showing how the number of polygons and vertices changes with different values of simplify.

[11]:
df = pd.DataFrame([polygon_metrics(m) for m in multipoly_list])
df = df.rename_axis("simplify")
df
[11]:
n_exterior_rings n_interior_rings n_exterior_vertices n_interior_vertices area
simplify
0 213 0 6621 0 1517.67225
1 160 0 6322 0 1514.93975
2 143 0 4795 0 1505.56895
3 127 0 4389 0 1498.73910
4 118 0 3403 0 1493.29985
5 109 0 2282 0 1489.30745
6 104 0 1641 0 1486.48475
7 100 0 1299 0 1481.45510
8 94 0 1073 0 1454.24045
9 90 0 911 0 1440.91335
10 84 0 652 0 1414.76040

Visualize#

Visualize how a few of the larger polygons vary with simplify.

[12]:
p0 = shapely.Point(-17, 50)

fig, ax = plt.subplots()
for i, multipoly in enumerate(multipoly_list):
    for p in multipoly.geoms:
        if p.contains(p0):
            ax.plot(*p.exterior.xy, label=f"simplify={i}")

ax.legend();
../_images/notebooks_polygons_19_0.png
[13]:
p0 = shapely.Point(-120, -54)

fig, ax = plt.subplots()
for i, multipoly in enumerate(multipoly_list):
    for p in multipoly.geoms:
        if p.contains(p0):
            ax.plot(*p.exterior.xy, label=f"simplify={i}")

ax.legend();
../_images/notebooks_polygons_20_0.png
[14]:
p0 = shapely.Point(55, 67)

fig, ax = plt.subplots()
for i, multipoly in enumerate(multipoly_list):
    for p in multipoly.geoms:
        if p.contains(p0):
            ax.plot(*p.exterior.xy, label=f"simplify={i}")

ax.legend();
../_images/notebooks_polygons_21_0.png

Experimental: convex hulls#

Use the boolean convex_hull query parameter to simplify polygons by computing their convex hulls. This parameter is experimental (false by default). When this parameter is enabled, the interiors parameter is automatically set to false.

This purpose of this parameter is to: - Further reduce the complexity of the polygon geometry. - Fill in “fjords” and other gaps in the polygon geometry. If polygons are regions to be avoided, points outside of a polygon but inside its convex hull may still be inaccessible for trajectory planning. By explicitly taking a convex hull, the flight planner may have an easier time creating optimal trajectories.

As seen in the table below, the number of polygons (n_exterior_rings) and the number of vertices are reduced when using the convex_hull parameter. On the other hand, the total footprint of the polygons is increased (as measured by the area metric).

[15]:
multipoly_list_ch = [get_shapely(convex_hull=True, simplify=s) for s in range(0, 11)]
[16]:
df1 = (
    df.drop(columns=["n_interior_rings", "n_interior_vertices"])
    .assign(convex_hull=False)
    .rename_axis("simplify", axis=0)
    .set_index("convex_hull", append=True)
)

df2 = pd.DataFrame([polygon_metrics(m) for m in multipoly_list_ch])
df2 = (
    df2.drop(columns=["n_interior_rings", "n_interior_vertices"])
    .assign(convex_hull=True)
    .rename_axis("simplify", axis=0)
    .set_index("convex_hull", append=True)
)

df = pd.concat([df1, df2], axis=0)
df.sort_index()
[16]:
n_exterior_rings n_exterior_vertices area
simplify convex_hull
0 False 213 6621 1517.67225
True 202 2187 2992.56350
1 False 160 6322 1514.93975
True 153 1870 2989.51695
2 False 143 4795 1505.56895
True 137 1602 2985.52305
3 False 127 4389 1498.73910
True 121 1341 2977.64030
4 False 118 3403 1493.29985
True 113 1102 2967.68565
5 False 109 2282 1489.30745
True 104 905 2951.80240
6 False 104 1641 1486.48475
True 99 809 2936.25635
7 False 100 1299 1481.45510
True 95 714 2914.85460
8 False 94 1073 1454.24045
True 89 622 2888.30430
9 False 90 911 1440.91335
True 86 556 2849.55900
10 False 84 652 1414.76040
True 82 428 2740.72190
[17]:
fig, ax = plt.subplots()

for poly in multipoly_list[6].geoms:
    ax.plot(*poly.exterior.xy, color="red", lw=1)

for poly in multipoly_list_ch[6].geoms:
    ax.plot(*poly.exterior.xy, color="purple", lw=1)

red_line = mlines.Line2D([], [], color="red", label="convex_hull = false")
purple_line = mlines.Line2D([], [], color="purple", label="convex_hull = true")
ax.legend(handles=[red_line, purple_line]);
../_images/notebooks_polygons_25_0.png
[ ]: