Inconsistent hole detection

[thabbott]

Description

Polygon generation with different thresholds can produce inconsistent sets of holes.

If a hole is present with a low threshold, it should also be present at a higher threshold provided the enclosing polygon remains intact. In some cases, however, holes appear only at intermediate thresholds.

In the example below, holes should be present in the large region of high EF at all thresholds, but appear only with thresholds near 50,000,000 J/m.

Details

• Version: v0.50.0
• Module: core.met
• OS: macOS

Steps to Reproduce

Download https://storage.googleapis.com/contrails-301217-public-data/sample.nc and run

import xarray as xr
import matplotlib.colors as colors
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import geojson
import numpy as np

from pycontrails.core import MetDataArray

def build_polygons(field, threshold):
   params = dict(
       fill_value=0.0,
       iso_value=threshold,
       min_area=0.3,
       epsilon=0.05,
       precision=2,
       interiors=True,
       convex_hull=False,
       include_altitude=True,
   )
   poly = field.to_polygon_feature(**params)
   return geojson.FeatureCollection(poly)

ds = xr.open_dataset("sample.nc")

plt.figure(figsize=(10, 3), dpi=300)
ax = plt.subplot(111, projection=ccrs.PlateCarree())
ax.coastlines(lw=0.1, color="gray")
im = ax.pcolormesh(ds["longitude"], ds["latitude"], ds["ef_per_m"].squeeze().T,
    shading="nearest", cmap="RdBu_r", norm=colors.SymLogNorm(linthresh=1e7))
plt.colorbar(im, label="EF (J m$^{-1}$)")

threshold = 3e7
for poly in build_polygons(
   MetDataArray(ds["ef_per_m"]),
   threshold
)["features"]["geometry"]["coordinates"]:
   for region in poly:
      v = np.array(region)
      ax.plot(v[:,0], v[:,1], transform=ccrs.Geodetic(), lw=0.5, color="black")

ax.set_title(f"{threshold} J/m")

plt.figure(figsize=(10, 3), dpi=300)
ax = plt.subplot(111, projection=ccrs.PlateCarree())
ax.coastlines(lw=0.1, color="gray")
im = ax.pcolormesh(ds["longitude"], ds["latitude"], ds["ef_per_m"].squeeze().T,
    shading="nearest", cmap="RdBu_r", norm=colors.SymLogNorm(linthresh=1e7))
plt.colorbar(im, label="EF (J m$^{-1}$)")

threshold = 5e7
for poly in build_polygons(
   MetDataArray(ds["ef_per_m"]),
   threshold
)["features"]["geometry"]["coordinates"]:
   for region in poly:
      v = np.array(region)
      ax.plot(v[:,0], v[:,1], transform=ccrs.Geodetic(), lw=0.5, color="black")

ax.set_title(f"{threshold} J/m")

plt.figure(figsize=(10, 3), dpi=300)
ax = plt.subplot(111, projection=ccrs.PlateCarree())
ax.coastlines(lw=0.1, color="gray")
im = ax.pcolormesh(ds["longitude"], ds["latitude"], ds["ef_per_m"].squeeze().T,
    shading="nearest", cmap="RdBu_r", norm=colors.SymLogNorm(linthresh=1e7))
plt.colorbar(im, label="EF (J m$^{-1}$)")

threshold = 7e7
for poly in build_polygons(
   MetDataArray(ds["ef_per_m"]),
   threshold
)["features"]["geometry"]["coordinates"]:
   for region in poly:
      v = np.array(region)
      ax.plot(v[:,0], v[:,1], transform=ccrs.Geodetic(), lw=0.5, color="black")

ax.set_title(f"{threshold} J/m")