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.

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

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

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

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.

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)

# 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

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

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}

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.

multipoly_list = [get_shapely(interiors=False, simplify=s) for s in range(0, 11)]

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

Metrics

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

df = pd.DataFrame([polygon_metrics(m) for m in multipoly_list])
df = df.rename_axis("simplify")
df

	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.

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();

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();

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();

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).

multipoly_list_ch = [get_shapely(convex_hull=True, simplify=s) for s in range(0, 11)]

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()

		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

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]);