Correct Ground Track Wrapping over Poles

[ptgrogan]

When wrapping polygons that cross the prime meridian over the north/south pole, make sure to split the polygon along the prime meridian to avoid the discontinuity along the anti-meridian after wrapping.

[ptgrogan]

Some more details on this issue -- TAT-C uses the sub-satellite point to determine the appropriate coordinate reference system (CRS) for the coverage buffering operation. However, the spatial extent of a coverage region for large-swath instruments extends far beyond the valid region for the universal transverse Mercator (UTM) zone when the sub-satellite point is close to a universal polar stenographic (UPS) zone but still within a UTM zone, causing highly distorted (wrapped on itself) polygons.

Here is an example visualization and supporting code:

from tatc import utils
from tatc.schemas import Instrument, Satellite, TwoLineElements

import geopandas as gpd

hirs = Instrument(
    name="HIRS",
    field_of_regard=utils.swath_width_to_field_of_regard(870000, 1820000), #1820 km swath at 870 km altitude
)
noaa_19 = Satellite(
    name="NOAA 19",
    orbit=TwoLineElements(
        tle=[
            "1 33591U 09005A   23230.15011991  .00000219  00000+0  14239-3 0  9995",
            "2 33591  99.0953 277.0538 0014800  70.7316 289.5454 14.12805463748692",
        ]
    ),
    instruments=[hirs],
)

# download ne_110m_admin_0_countries.zip from https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip
worldmap = gpd.read_file("ne_110m_admin_0_countries.zip")

from tatc.analysis import collect_ground_track
from datetime import datetime, timedelta, timezone
import pandas as pd

start = datetime(year=2023, month=8, day=10, hour=0, minute=30, tzinfo=timezone.utc)
duration = timedelta(minutes=10)
time_step = timedelta(seconds=10)

dates = pd.date_range(
    start, 
    start + duration,
    freq=time_step
)

# cylindrical crs
results_4087 = collect_ground_track(noaa_19, hirs, dates, crs="EPSG:4087")
results_4087["crs"] = "EPSG:4087"

# web mercator crs
results_3857 = collect_ground_track(noaa_19, hirs, dates, crs="EPSG:3857")
results_3857["crs"] = "EPSG:3857"

# universal polar sterographic crs
results_32661 = collect_ground_track(noaa_19, hirs, dates, crs="EPSG:5041")
results_32661["crs"] = "EPSG:5041"

# universal transverse mercator crs (uses universal polar sterographic crs for polar sub-satellite points)
results_utm = collect_ground_track(noaa_19, hirs, dates, crs="utm")
results_utm["crs"] = "utm"

results = pd.concat([results_4087, results_3857, results_32661, results_utm])

from tatc.analysis import collect_orbit_track

# sub-satellite points natively use WGS84 crs
results_ssp = collect_orbit_track(noaa_19, hirs, dates)
results_ssp["crs"] = "EPSG:4326"

# extract the crs used if the user request the utm crs
from tatc.analysis.track import _get_utm_epsg_code
results_ssp["utm"] = results_ssp.apply(lambda row: _get_utm_epsg_code(row.geometry), axis=1)

import matplotlib.animation as animation
import matplotlib.pyplot as plt
import geopandas as gpd

def animate_frame(ax, crs, frame):
    ax.clear()
    filtered_results = results[
        (results["time"] >= start + frame*time_step) 
        & (results["time"] < start + (frame + 1)*time_step)
        & (results["crs"] == crs)
    ]
    filtered_results.plot(ax=ax)
    worldmap.dissolve().boundary.plot(ax=ax, color="black", lw=0.5)
    if crs == "utm":
        ax.set_title(f"utm ({results_ssp[results_ssp.time==filtered_results.iloc[0].time].iloc[0].utm})")
    else:
        ax.set_title(crs)
    ax.set_xlim((-180,180))
    ax.set_ylim((-90,90))
    ax.set_aspect('equal')

fig, axs = plt.subplots(2,2, figsize=(8,6))

def animate(frame):
    animate_frame(axs[0,0], "EPSG:4087", frame)
    animate_frame(axs[0,1], "EPSG:3857", frame)
    animate_frame(axs[1,0], "EPSG:5041", frame)
    animate_frame(axs[1,1], "utm", frame)
    plt.suptitle(f"Frame {frame}")

ani = animation.FuncAnimation(fig, animate, frames=int(duration/time_step), interval=200, blit=False)
ani.save("tatc-error.gif", dpi=300, writer=animation.ImageMagickWriter(fps=5))