example: select obs based on a netcdf file, e.g. mpas region

[hkershaw-brown]

This might be a nice example to have.
Read in mpas file, create polygon of boundary, trim obs_seq

Playing around (note this is a super bad example since I have an ocean obs_seq).

import xarray as xr
from shapely.geometry import MultiPoint, Polygon
import numpy as np

# Open the MPAS restart file
ds = xr.open_dataset('../data/wofs_mpas.afterda.nc')

# Get all vertex coordinates
lon_vertices = ds.lonVertex.values
lat_vertices = ds.latVertex.values

# Build list of (lon, lat) tuples
all_vertex_coords = list(zip(lon_vertices, lat_vertices))

# Create a MultiPoint object and get its convex hull
boundary_polygon = MultiPoint(all_vertex_coords).convex_hull

import matplotlib.pyplot as plt
import numpy as np

# Convert boundary_polygon coordinates from radians to degrees
if boundary_polygon.geom_type == 'Polygon':
    x_rad, y_rad = boundary_polygon.exterior.xy
    x_deg = np.degrees(x_rad)
    y_deg = np.degrees(y_rad)
else:
    x_deg, y_deg = [], []
    for geom in boundary_polygon.geoms:
        x_rad, y_rad = geom.exterior.xy
        x_deg.extend(np.degrees(x_rad))
        y_deg.extend(np.degrees(y_rad))

fig, ax = plt.subplots(figsize=(10, 8))

# Plot all observation points from mom6.df (GeoDataFrame)
gdf.plot(ax=ax, color='gray', alpha=0.5, markersize=5, label='All Observations')

# Plot the MPAS grid boundary polygon in degrees
ax.plot(x_deg, y_deg, color='green', linewidth=2, label='MPAS Grid Boundary (deg)')

ax.set_xlabel('Longitude (deg)')
ax.set_ylabel('Latitude (deg)')
ax.set_title('All Observations and MPAS Grid Boundary (Degrees)')
ax.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

I thnk geopandas + shapely will take care of the masking, geometry.within

[hkershaw-brown]

see NCAR/DART#935