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 = "https://api.contrails.org/v0"
api_key = os.environ["CONTRAILS_API_KEY"]
headers = {"x-api-key": api_key}
# Confirm credentials are valid
r = requests.get(f"{URL}/auth/validate-key", headers=headers)
print(f"HTTP Response Code: {r.status_code} {r.reason}")
HTTP Response Code: 200 OK
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
: Setformat=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
: Setinteriors=true
to return interior polygons.simplify
: Used to control the extent to which polygon geometry is simplified.convex_hull
: Setconvex_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)

[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}

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)

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

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

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

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

[ ]: