From 58e2bdc402f6f2eb6d1f523066b1b4ed0ebf091f Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 25 Jun 2026 12:58:57 -0800 Subject: [PATCH 01/17] Clean up structure of AIS reader.py file + transition towards accepting paths --- tests/utils/ais/test_alignment.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/tests/utils/ais/test_alignment.py b/tests/utils/ais/test_alignment.py index dbcea43..11e70b1 100644 --- a/tests/utils/ais/test_alignment.py +++ b/tests/utils/ais/test_alignment.py @@ -19,8 +19,15 @@ from nps_active_space.utils.time_utils import site_timezone_name, utc_naive_to_site_naive REPO = Path(__file__).resolve().parents[3] -AIS_MAY24 = REPO / "example_data" / "AIS" / "MXAK-AIS-GLBA-20240524.csv" -NVSPL_DIR = REPO / "example_data" / "NVSPL" / "GLBALSTL_2024" +AIS_MAY24 = REPO / "example_data" / "MXAK-AIS-GLBA" / "MXAK-AIS-GLBA-20240524.csv" +NVSPL_DIR = ( + REPO + / "example_data" + / "nvspl_archive" + / "2024 GLBALSTL Lower South Tidal Inlet" + / "01 DATA" + / "NVSPL" +) SHIP_VISITS = REPO / "example_data" / "GLBALSTL_ship_visits.csv" GLBALSTL_LAT, GLBALSTL_LON = 58.78, -136.32 @@ -29,7 +36,7 @@ # Same transit scale as MXAK event splitting; used for CPA windows and loose AIS timing. CPA_WINDOW_MINUTES = EVENT_GAP_SECONDS // 60 AIS_CPA_TOLERANCE_SEC = 2 * EVENT_GAP_SECONDS -UTC_MISALIGNMENT_FLOOR_SEC = 60 * 60 +UTC_MISALIGNMENT_FLOOR_SEC = 3600 NVSPL_PEAK_TOLERANCE_DB = 5 MAY24_MMSIS = (367365630, 246648000, 367578110) From f9e25cca3a3bf2871bf746aa64d625c51e332864 Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 11:48:50 -0800 Subject: [PATCH 02/17] test vis changes for AIS tracks --- nps_active_space/scripts/viz.py | 185 ++++++++++++++++++++++++++++++-- tests/scripts/test_viz.py | 79 ++++++++++++++ 2 files changed, 256 insertions(+), 8 deletions(-) create mode 100644 tests/scripts/test_viz.py diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index 6abc109..8d526b8 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -8,9 +8,10 @@ import pyproj import sys from shapely.geometry import box, Polygon, MultiPolygon, LineString, MultiLineString -from nps_active_space.utils.helpers import get_deployment, load_annotations, load_DEM, load_layered_activespace, load_studyarea +from nps_active_space.utils.ais import query_ais_mxak +from nps_active_space.utils.helpers import get_deployment, get_elevation, load_annotations, load_DEM, load_layered_activespace, load_studyarea from nps_active_space.scripts.run_audible_transits import AudibleTransits -from nps_active_space.utils.models import Annotations +from nps_active_space.utils.models import Annotations, Tracks from nps_active_space.utils.computation import NMSIM_bbox_utm import nps_active_space.utils.config as cfg import argparse @@ -43,12 +44,50 @@ def active_to_linestrings(active): return linestrings +def densify_linestring(linestring: LineString, step_m: float) -> LineString: + """Insert vertices along a line so elevation can be sampled between sparse AIS points.""" + if linestring.is_empty or linestring.length <= step_m: + return linestring + distances = np.arange(0, linestring.length, step_m) + if distances[-1] < linestring.length: + distances = np.append(distances, linestring.length) + return LineString([linestring.interpolate(float(d)) for d in distances]) + + +def sea_surface_z_profile( + linestring: LineString, + dem, + crs: str, + *, + offset_m: float = 5.0, + densify_step_m: float = 100.0, +) -> tuple[LineString, np.ndarray]: + """Return a densified line and z values hugging the DEM water surface (≈0 m).""" + line = densify_linestring(linestring, densify_step_m) + to_wgs84 = pyproj.Transformer.from_crs(crs, "epsg:4326", always_xy=True) + z_vals = [] + for x, y in np.array(line.coords)[:, :2]: + lon, lat = to_wgs84.transform(x, y) + dem_z = float(get_elevation(dem, lon, lat)) + # NMSIM DEM uses 0 m over water; keep tracks slightly above the mesh to avoid z-fighting. + z_vals.append(max(dem_z, 0.0) + offset_m) + return line, np.asarray(z_vals) + + def create_polyline_3d(linestring, z=None): coords = np.array(linestring.coords) + xy = coords[:, :2] if z is not None: - coords = np.column_stack((coords, np.full(coords.shape[0], z))) + z_arr = np.atleast_1d(np.asarray(z, dtype=float)) + if z_arr.size != coords.shape[0]: + raise ValueError( + f"z must be scalar or have one value per vertex ({coords.shape[0]}), got {z_arr.size}" + ) + coords = np.column_stack((xy, z_arr)) + elif coords.shape[1] == 2: + coords = np.column_stack((xy, np.zeros(coords.shape[0]))) assert coords.shape[1] == 3 - coords[:,2] = np.clip(coords[:,2], 0, 10000) # reduce magnitude of erroneous elevations + coords[:, 2] = np.clip(coords[:, 2], 0, 10000) # reduce magnitude of erroneous elevations n_points = coords.shape[0] lines = np.hstack([[n_points], np.arange(n_points)]) poly = pv.PolyData(coords) @@ -56,6 +95,14 @@ def create_polyline_3d(linestring, z=None): return poly +def track_points_to_linestring(track_points: gpd.GeoDataFrame) -> LineString: + """Connect ordered track points into a 2D line in the plot CRS.""" + coords = [(geom.x, geom.y) for geom in track_points.geometry] + if len(coords) < 2: + return LineString(coords) + return LineString(coords) + + # main class ===================================================== @@ -67,11 +114,15 @@ class Visualizer(): audible_annotation_color = "deepskyblue" inaudible_annotation_color = "red" audible_transits_color = "purple" + vessel_track_color = "cyan" z_scale_toggle_color = "black" # button toggling z scale + sea_surface_offset_m = 5.0 + sea_surface_densify_step_m = 100.0 def __init__(self, unit, site, year, env, do_active=False, gain=None, - do_annots=False, do_transits=False, + do_annots=False, do_transits=False, do_vessels=False, annotation_file=None, audible_transits_pkl=None, + vessel_start_date=None, vessel_end_date=None, terraced=False, fill_layers=False, max_tracks=1000, ): # class metadata @@ -99,6 +150,8 @@ def __init__(self, unit, site, year, env, do_active=False, gain=None, self.plot_annotations(annotation_file) if do_transits: self.plot_audible_transits(audible_transits_pkl) + if do_vessels: + self.plot_vessel_tracks(vessel_start_date, vessel_end_date) # configure plot parameters and display self.plotter.camera_position = "xz" @@ -111,7 +164,8 @@ def __init__(self, unit, site, year, env, do_active=False, gain=None, def plot_dem(self, show_scalar_bar=False): # load DEM - dem = load_DEM(self.project_dir, self.unit, self.site) + dem = load_DEM(self.project_dir, self.unit, self.site) + self.dem = dem data = dem.read(1) if dem.nodata is not None: data[data == dem.nodata] = 0 @@ -266,6 +320,44 @@ def toggle(flag): color_on=self.activespace_color ) + def _annotation_z_values(self, linestring, offset_m: float = 2.0) -> np.ndarray: + """Sample DEM elevation when track z is at or below the surface (e.g. vessel annotations).""" + coords = np.array(linestring.coords) + vertex_z = [coord[2] if len(coord) > 2 else 0.0 for coord in coords] + if all(z <= 0 for z in vertex_z): + _, z_vals = sea_surface_z_profile( + linestring, + self.dem, + self.crs, + offset_m=offset_m, + densify_step_m=self.sea_surface_densify_step_m, + ) + return z_vals + + to_wgs84 = pyproj.Transformer.from_crs(self.crs, "epsg:4326", always_xy=True) + z_vals = [] + for coord in coords: + x, y = coord[0], coord[1] + z = coord[2] if len(coord) > 2 else 0.0 + lon, lat = to_wgs84.transform(x, y) + dem_z = float(get_elevation(self.dem, lon, lat)) + if z <= 0 or z < dem_z: + z_vals.append(max(dem_z, 0.0) + offset_m) + else: + z_vals.append(z) + return np.asarray(z_vals) + + def _sea_surface_polyline(self, linestring: LineString) -> pv.PolyData: + """Build a 3D polyline that follows the local water/ground surface.""" + line, z_vals = sea_surface_z_profile( + linestring, + self.dem, + self.crs, + offset_m=self.sea_surface_offset_m, + densify_step_m=self.sea_surface_densify_step_m, + ) + return create_polyline_3d(line, z=z_vals) + def plot_annotations(self, annotation_file=None): # load annotations if annotation_file is None: @@ -296,7 +388,14 @@ def plot_annotations(self, annotation_file=None): inaudible_actors = [] for _, annot in annotations.iterrows(): color = "deepskyblue" if annot["audible"] else "red" - polyline = create_polyline_3d(annot["geometry"]) + geometry = annot["geometry"] + coords = np.array(geometry.coords) + vertex_z = [c[2] if len(c) > 2 else 0.0 for c in coords] + if all(z <= 0 for z in vertex_z): + polyline = self._sea_surface_polyline(geometry) + else: + z_vals = self._annotation_z_values(geometry) + polyline = create_polyline_3d(geometry, z=z_vals) actor = self.plotter.add_mesh(polyline, color=color, point_size=2, line_width=2) if annot["audible"]: audible_actors.append(actor) @@ -373,6 +472,70 @@ def toggle(flag): color_on="purple" ) + def plot_vessel_tracks(self, start_date: str | None = None, end_date: str | None = None): + """Plot MXAK AIS vessel transits at the local sea surface.""" + start_date = start_date or f"{self.year}-01-01" + end_date = end_date or f"{self.year}-12-31" + + try: + ais_path = cfg.read("data", "ais") + except KeyError: + print("No [data] ais path in config, skipping vessel tracks.") + return + + print(f"Querying AIS tracks from {start_date} to {end_date}") + try: + raw_tracks = query_ais_mxak( + ais_path=ais_path, + start_date=start_date, + end_date=end_date, + mask=self.study_area, + ) + except AssertionError as exc: + print(f"No AIS tracks loaded: {exc}") + return + tracks = Tracks(raw_tracks, id_col="event_id", datetime_col="TIME", z_col="altitude") + tracks = tracks.to_crs(self.crs) + + track_ids = tracks["track_id"].drop_duplicates() + print(f"{len(track_ids)} vessel transits ({len(tracks)} points)") + if len(track_ids) > self.max_tracks: + print(f"More than {self.max_tracks}, sampling") + selected = track_ids.sample(self.max_tracks, replace=False, random_state=3) + tracks = tracks[tracks["track_id"].isin(selected)] + print(f"Showing {selected.nunique()} transits") + + actors = [] + for track_id, group in tracks.groupby("track_id", sort=False): + group = group.sort_values("point_dt") + line = track_points_to_linestring(group) + if line.is_empty or line.length == 0: + continue + polyline = self._sea_surface_polyline(line) + actor = self.plotter.add_mesh( + polyline, + color=self.vessel_track_color, + point_size=2, + line_width=3, + ) + actors.append(actor) + + if not actors: + print("No vessel tracks to plot.") + return + + def toggle(flag): + for actor in actors: + actor.SetVisibility(flag) + + self.plotter.add_checkbox_button_widget( + callback=toggle, + value=True, + position=(10, 220), + size=25, + color_on=self.vessel_track_color, + ) + def setup_z_scale(self): self.plotter.set_scale(1, 1, 2) # easier to make out elevation differences @@ -408,6 +571,8 @@ def toggle_z_scale(flag): help="If included, load and plot annotations") parser.add_argument("-t", "--audible-transits", action="store_true", help="If included, load and plot audible transits") + parser.add_argument("-v", "--vessels", action="store_true", + help="If included, load and plot MXAK AIS vessel tracks at sea level") parser.add_argument("--all", action="store_true", help="Load everything, shorthand for --active-space --annotations --audible-transits") @@ -416,6 +581,8 @@ def toggle_z_scale(flag): help="Maximum number of annotation tracks or audible transits to show.") parser.add_argument("--annotation-file", help="Path to .geojson file to load annotations from, if not the default.") parser.add_argument("--transits-pkl", help="Path to .pkl file to load audible transits from, if not the default.") + parser.add_argument("--start-date", help="AIS query start date (YYYY-MM-DD). Default: Jan 1 of deployment year.") + parser.add_argument("--end-date", help="AIS query end date (YYYY-MM-DD). Default: Dec 31 of deployment year.") parser.add_argument("--terraced", action="store_true", help="If included, render the active space as the terraced surface instead of contours.") parser.add_argument("--fill-layers", action="store_true", @@ -430,9 +597,11 @@ def toggle_z_scale(flag): do_active = args.active_space or args.all do_annotations = args.annotations or args.all do_transits = args.audible_transits or args.all + do_vessels = args.vessels Visualizer(unit, site, year, args.environment, do_active, args.gain, - do_annotations, do_transits, + do_annotations, do_transits, do_vessels, args.annotation_file, args.transits_pkl, + args.start_date, args.end_date, args.terraced, args.fill_layers, args.max_tracks) \ No newline at end of file diff --git a/tests/scripts/test_viz.py b/tests/scripts/test_viz.py new file mode 100644 index 0000000..3e7b68f --- /dev/null +++ b/tests/scripts/test_viz.py @@ -0,0 +1,79 @@ +import numpy as np +import pytest +from shapely.geometry import LineString, Point + +from nps_active_space.scripts.viz import ( + create_polyline_3d, + densify_linestring, + sea_surface_z_profile, + track_points_to_linestring, +) + + +class TestDensifyLinestring: + def test_short_line_unchanged(self): + line = LineString([(0, 0), (50, 0)]) + assert densify_linestring(line, step_m=100) == line + + def test_long_line_gets_intermediate_vertices(self): + line = LineString([(0, 0), (250, 0)]) + dense = densify_linestring(line, step_m=100) + assert len(dense.coords) == 4 + assert dense.coords[1] == pytest.approx((100.0, 0.0)) + + +class TestSeaSurfaceZProfile: + def test_water_pixels_use_offset_above_zero(self, monkeypatch): + line = LineString([(0, 0), (250, 0)]) + + def fake_elevation(_dem, _lon, _lat): + return 0.0 + + monkeypatch.setattr( + "nps_active_space.scripts.viz.get_elevation", + fake_elevation, + ) + dense, z_vals = sea_surface_z_profile( + line, + dem=object(), + crs="epsg:32608", + offset_m=5.0, + densify_step_m=100.0, + ) + assert len(dense.coords) == 4 + assert z_vals.tolist() == [5.0, 5.0, 5.0, 5.0] + + def test_negative_dem_clamped_to_sea_level(self, monkeypatch): + line = LineString([(0, 0), (1, 0)]) + + monkeypatch.setattr( + "nps_active_space.scripts.viz.get_elevation", + lambda _dem, _lon, _lat: -12.0, + ) + _, z_vals = sea_surface_z_profile( + line, + dem=object(), + crs="epsg:32608", + offset_m=2.0, + densify_step_m=100.0, + ) + assert z_vals.tolist() == [2.0, 2.0] + + +class TestCreatePolyline3d: + def test_rejects_mismatched_z_length(self): + line = LineString([(0, 0), (1, 1), (2, 2)]) + with pytest.raises(ValueError, match="one value per vertex"): + create_polyline_3d(line, z=np.array([1.0, 2.0])) + + +class TestTrackPointsToLinestring: + def test_builds_line_from_points(self): + import geopandas as gpd + + gdf = gpd.GeoDataFrame( + geometry=[Point(0, 0), Point(1, 0), Point(2, 0)], + crs="epsg:32608", + ) + line = track_points_to_linestring(gdf) + assert len(line.coords) == 3 From fe99666e0a53b26c1a9472bda49e7fa86732ca6c Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 11:52:46 -0800 Subject: [PATCH 03/17] Type hints and better CLI arg validation --- nps_active_space/scripts/viz.py | 180 +++++++++++++++++++++++++------- tests/scripts/test_viz.py | 63 +++++++++++ 2 files changed, 205 insertions(+), 38 deletions(-) diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index 8d526b8..f0526aa 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -1,24 +1,84 @@ +from __future__ import annotations + +import argparse +import glob import os +import sys +from datetime import datetime +from pathlib import Path + +import geopandas as gpd import numpy as np import pandas as pd -import geopandas as gpd -import glob +import pyproj import pyvista as pv import rasterio -import pyproj -import sys -from shapely.geometry import box, Polygon, MultiPolygon, LineString, MultiLineString +from shapely.geometry import LineString, MultiLineString, MultiPolygon, Polygon, box from nps_active_space.utils.ais import query_ais_mxak from nps_active_space.utils.helpers import get_deployment, get_elevation, load_annotations, load_DEM, load_layered_activespace, load_studyarea from nps_active_space.scripts.run_audible_transits import AudibleTransits from nps_active_space.utils.models import Annotations, Tracks from nps_active_space.utils.computation import NMSIM_bbox_utm import nps_active_space.utils.config as cfg -import argparse + +# CLI argument validation ===================================== + +def parse_deployment(name: str) -> tuple[str, str, int]: + """Parse deployment name like DENATRLA2024 into unit, site, and year.""" + if len(name) < 9: + raise argparse.ArgumentTypeError( + f"deployment must be at least 9 characters " + f"(4-char unit + site code + 4-digit year), got {len(name)}: {name!r}" + ) + unit, site, year_str = name[:4], name[4:-4], name[-4:] + if not site: + raise argparse.ArgumentTypeError( + f"deployment must include a site code between unit and year, got {name!r}" + ) + if len(year_str) != 4 or not year_str.isdigit(): + raise argparse.ArgumentTypeError( + f"deployment year must be 4 digits, got {year_str!r} in {name!r}" + ) + return unit, site, int(year_str) + + +def parse_iso_date(value: str, *, arg_name: str) -> str: + """Validate YYYY-MM-DD date string.""" + try: + datetime.strptime(value, "%Y-%m-%d") + except ValueError: + raise argparse.ArgumentTypeError( + f"{arg_name}: expected YYYY-MM-DD date, got {value!r}" + ) from None + return value + + +def parse_existing_file(path: str, *, arg_name: str) -> str: + """Validate that a CLI path refers to an existing file.""" + file_path = Path(path) + if not file_path.is_file(): + raise argparse.ArgumentTypeError(f"{arg_name}: file not found: {path}") + return path + + +def parse_max_tracks(value: str) -> int: + """Validate positive integer for --max-tracks.""" + try: + n = int(value) + except ValueError: + raise argparse.ArgumentTypeError( + f"--max-tracks must be a positive integer, got {value!r}" + ) from None + if n <= 0: + raise argparse.ArgumentTypeError( + f"--max-tracks must be a positive integer, got {n}" + ) + return n + # helper functions ============================================ -def polygon_to_mesh(polygon, z): +def polygon_to_mesh(polygon: Polygon, z: float) -> pv.PolyData: exterior = np.array(polygon.exterior.coords) points = np.c_[exterior[:, 0], exterior[:, 1], np.full(len(exterior), z)] poly = pv.PolyData(points) @@ -26,7 +86,7 @@ def polygon_to_mesh(polygon, z): return poly.triangulate() -def active_to_polys(active): +def active_to_polys(active: gpd.GeoDataFrame) -> list[Polygon]: geometry = active.geometry.iloc[0] if isinstance(geometry, Polygon): polys = [geometry] @@ -35,7 +95,7 @@ def active_to_polys(active): return polys -def active_to_linestrings(active): +def active_to_linestrings(active: gpd.GeoDataFrame) -> list[LineString]: geometry = active.boundary.iloc[0] if isinstance(geometry, MultiLineString): linestrings = geometry.geoms @@ -74,7 +134,10 @@ def sea_surface_z_profile( return line, np.asarray(z_vals) -def create_polyline_3d(linestring, z=None): +def create_polyline_3d( + linestring: LineString, + z: float | np.ndarray | None = None, +) -> pv.PolyData: coords = np.array(linestring.coords) xy = coords[:, :2] if z is not None: @@ -119,12 +182,25 @@ class Visualizer(): sea_surface_offset_m = 5.0 sea_surface_densify_step_m = 100.0 - def __init__(self, unit, site, year, env, do_active=False, gain=None, - do_annots=False, do_transits=False, do_vessels=False, - annotation_file=None, audible_transits_pkl=None, - vessel_start_date=None, vessel_end_date=None, - terraced=False, fill_layers=False, max_tracks=1000, - ): + def __init__( + self, + unit: str, + site: str, + year: int, + env: str, + do_active: bool = False, + gain: float | None = None, + do_annots: bool = False, + do_transits: bool = False, + do_vessels: bool = False, + annotation_file: str | None = None, + audible_transits_pkl: str | None = None, + vessel_start_date: str | None = None, + vessel_end_date: str | None = None, + terraced: bool = False, + fill_layers: bool = False, + max_tracks: int = 1000, + ) -> None: # class metadata self.unit = unit self.site = site @@ -162,7 +238,7 @@ def __init__(self, unit, site, year, env, do_active=False, gain=None, self.plotter.show_axes() self.plotter.show() - def plot_dem(self, show_scalar_bar=False): + def plot_dem(self, show_scalar_bar: bool = False) -> None: # load DEM dem = load_DEM(self.project_dir, self.unit, self.site) self.dem = dem @@ -192,23 +268,23 @@ def plot_dem(self, show_scalar_bar=False): # add mesh to plot self.plotter.add_mesh(mesh, scalars="elevation", cmap="gist_earth", show_scalar_bar=show_scalar_bar) - def plot_point(self, x, y, z, color="white"): + def plot_point(self, x: float, y: float, z: float, color: str = "white") -> None: point = pv.PolyData(np.array([[x, y, z]])) self.plotter.add_mesh(point, color=color, point_size=10, render_points_as_spheres=True) - def plot_mic(self): + def plot_mic(self) -> None: mic = get_deployment(self.project_dir, self.unit, self.site, self.year) mic = mic.to_crs(self.crs) self.plot_point(mic.x, mic.y, mic.z, self.mic_color) - def plot_activespace(self, terraced=False, gain=None): + def plot_activespace(self, terraced: bool = False, gain: float | None = None) -> None: csv_3d_fits = f"{self.project_dir}/fits.csv" fit_results = pd.read_csv(csv_3d_fits, index_col="Designator") if gain is not None: print(f"Using gain {gain}dB") elif self.deployment in fit_results.index: - gain = float(fit_results.loc[unit+site+year, "1/3rd Octave Gain (F1)"]) + gain = float(fit_results.loc[self.deployment, "1/3rd Octave Gain (F1)"]) else: print(f"No fitted active space gain found in {csv_3d_fits}, skipping active space.") return @@ -221,7 +297,7 @@ def plot_activespace(self, terraced=False, gain=None): else: self.plot_contoured_activespace(active_3d) - def plot_contoured_activespace(self, active_3d): + def plot_contoured_activespace(self, active_3d) -> None: layer_checkboxes = [] layer_callbacks = [] for i, (active_z, active) in enumerate(active_3d.activespaces.items()): @@ -245,7 +321,7 @@ def toggle_all_actives(flag): color_on=self.activespace_color ) - def plot_active_layer(self, active_layer, elevation, i=0): + def plot_active_layer(self, active_layer, elevation: float, i: int = 0): poly_actor = None if self.fill_layers: meshes = [] @@ -276,7 +352,9 @@ def toggle(flag): return checkbox, toggle - def plot_terraced_activespace(self, active_3d, layer_thickness=300, opacity=1): + def plot_terraced_activespace( + self, active_3d, layer_thickness: int = 300, opacity: float = 1 + ) -> None: meshes = [] layers = list(active_3d.activespaces.items()) for i in range(len(layers)): @@ -358,7 +436,7 @@ def _sea_surface_polyline(self, linestring: LineString) -> pv.PolyData: ) return create_polyline_3d(line, z=z_vals) - def plot_annotations(self, annotation_file=None): + def plot_annotations(self, annotation_file: str | None = None) -> None: # load annotations if annotation_file is None: annotations = load_annotations(self.project_dir, self.unit, self.site, self.year, only_valid=True) @@ -425,7 +503,7 @@ def toggle_inaudible(flag): color_on="red" ) - def plot_audible_transits(self, audible_transits_pkl=None): + def plot_audible_transits(self, audible_transits_pkl: str | None = None) -> None: if audible_transits_pkl: print(f"Loading audible transits from {audible_transits_pkl}") listener = AudibleTransits.from_pickle(audible_transits_pkl) @@ -433,7 +511,7 @@ def plot_audible_transits(self, audible_transits_pkl=None): print("Loading audible transits") matches = glob.glob(os.path.join( self.project_dir, self.unit+self.site, "Output_Data", "AUDIBLE_TRANSITS", - f"3D*{year}-01-01*Active Space {self.year}*", "AudibleTransits_object.pkl")) + f"3D*{self.year}-01-01*Active Space {self.year}*", "AudibleTransits_object.pkl")) if len(matches) == 0: print("No audible transits pkl file found found") return @@ -472,7 +550,9 @@ def toggle(flag): color_on="purple" ) - def plot_vessel_tracks(self, start_date: str | None = None, end_date: str | None = None): + def plot_vessel_tracks( + self, start_date: str | None = None, end_date: str | None = None + ) -> None: """Plot MXAK AIS vessel transits at the local sea surface.""" start_date = start_date or f"{self.year}-01-01" end_date = end_date or f"{self.year}-12-31" @@ -536,7 +616,7 @@ def toggle(flag): color_on=self.vessel_track_color, ) - def setup_z_scale(self): + def setup_z_scale(self) -> None: self.plotter.set_scale(1, 1, 2) # easier to make out elevation differences # allow toggling though @@ -559,7 +639,11 @@ def toggle_z_scale(flag): # required args # we pass deployment as a single positional arg instead of -u DENA -s TRLA -y 2024, # because it saves a lot of keystrokes, and this script is often used frequently - parser.add_argument("deployment", help="Deployment name, e.g. DENATRLA2024") + parser.add_argument( + "deployment", + type=parse_deployment, + help="Deployment name, e.g. DENATRLA2024", + ) parser.add_argument("-e", "--environment", default="DENA_streamline", help="Config environment name, e.g. DENA_streamline. Default is 'DENA_streamline'") @@ -577,12 +661,33 @@ def toggle_z_scale(flag): help="Load everything, shorthand for --active-space --annotations --audible-transits") # uncommon / special use case args - parser.add_argument("-m", "--max-tracks", default=500, - help="Maximum number of annotation tracks or audible transits to show.") - parser.add_argument("--annotation-file", help="Path to .geojson file to load annotations from, if not the default.") - parser.add_argument("--transits-pkl", help="Path to .pkl file to load audible transits from, if not the default.") - parser.add_argument("--start-date", help="AIS query start date (YYYY-MM-DD). Default: Jan 1 of deployment year.") - parser.add_argument("--end-date", help="AIS query end date (YYYY-MM-DD). Default: Dec 31 of deployment year.") + parser.add_argument( + "-m", + "--max-tracks", + type=parse_max_tracks, + default=500, + help="Maximum number of annotation tracks or audible transits to show.", + ) + parser.add_argument( + "--annotation-file", + type=lambda p: parse_existing_file(p, arg_name="--annotation-file"), + help="Path to .geojson file to load annotations from, if not the default.", + ) + parser.add_argument( + "--transits-pkl", + type=lambda p: parse_existing_file(p, arg_name="--transits-pkl"), + help="Path to .pkl file to load audible transits from, if not the default.", + ) + parser.add_argument( + "--start-date", + type=lambda d: parse_iso_date(d, arg_name="--start-date"), + help="AIS query start date (YYYY-MM-DD). Default: Jan 1 of deployment year.", + ) + parser.add_argument( + "--end-date", + type=lambda d: parse_iso_date(d, arg_name="--end-date"), + help="AIS query end date (YYYY-MM-DD). Default: Dec 31 of deployment year.", + ) parser.add_argument("--terraced", action="store_true", help="If included, render the active space as the terraced surface instead of contours.") parser.add_argument("--fill-layers", action="store_true", @@ -590,8 +695,7 @@ def toggle_z_scale(flag): args = parser.parse_args() - usy = args.deployment - unit, site, year = usy[:4], usy[4:-4], usy[-4:] + unit, site, year = args.deployment print(unit, site, year) do_active = args.active_space or args.all diff --git a/tests/scripts/test_viz.py b/tests/scripts/test_viz.py index 3e7b68f..1b2d637 100644 --- a/tests/scripts/test_viz.py +++ b/tests/scripts/test_viz.py @@ -1,3 +1,6 @@ +import argparse +from pathlib import Path + import numpy as np import pytest from shapely.geometry import LineString, Point @@ -5,11 +8,71 @@ from nps_active_space.scripts.viz import ( create_polyline_3d, densify_linestring, + parse_deployment, + parse_existing_file, + parse_iso_date, + parse_max_tracks, sea_surface_z_profile, track_points_to_linestring, ) +class TestParseDeployment: + def test_valid_deployment(self): + assert parse_deployment("DENATRLA2024") == ("DENA", "TRLA", 2024) + + def test_rejects_too_short(self): + with pytest.raises(argparse.ArgumentTypeError, match="at least 9 characters"): + parse_deployment("DENA2024") + + def test_rejects_unit_year_only(self): + with pytest.raises(argparse.ArgumentTypeError, match="at least 9 characters"): + parse_deployment("DENA2024") + + def test_rejects_non_digit_year(self): + with pytest.raises(argparse.ArgumentTypeError, match="4 digits"): + parse_deployment("DENATRLA20XX") + + +class TestParseIsoDate: + def test_valid_date(self): + assert parse_iso_date("2024-06-18", arg_name="--start-date") == "2024-06-18" + + def test_rejects_bad_format(self): + with pytest.raises(argparse.ArgumentTypeError, match="YYYY-MM-DD"): + parse_iso_date("06/18/2024", arg_name="--start-date") + + def test_includes_arg_name_in_error(self): + with pytest.raises(argparse.ArgumentTypeError, match="--end-date"): + parse_iso_date("not-a-date", arg_name="--end-date") + + +class TestParseExistingFile: + def test_existing_file(self, tmp_path: Path): + file_path = tmp_path / "annotations.geojson" + file_path.write_text("{}") + assert parse_existing_file(str(file_path), arg_name="--annotation-file") == str( + file_path + ) + + def test_missing_file(self, tmp_path: Path): + with pytest.raises(argparse.ArgumentTypeError, match="file not found"): + parse_existing_file(str(tmp_path / "missing.geojson"), arg_name="--transits-pkl") + + +class TestParseMaxTracks: + def test_valid_positive_int(self): + assert parse_max_tracks("500") == 500 + + def test_rejects_zero(self): + with pytest.raises(argparse.ArgumentTypeError, match="positive integer"): + parse_max_tracks("0") + + def test_rejects_non_integer(self): + with pytest.raises(argparse.ArgumentTypeError, match="positive integer"): + parse_max_tracks("abc") + + class TestDensifyLinestring: def test_short_line_unchanged(self): line = LineString([(0, 0), (50, 0)]) From 64da1ee09d40282dc33b9922377c59660633d31c Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 11:56:54 -0800 Subject: [PATCH 04/17] test cursor fix for not displaying AIS tracks --- nps_active_space/scripts/viz.py | 117 ++++++++++++++++++++++---------- 1 file changed, 83 insertions(+), 34 deletions(-) diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index f0526aa..70802b3 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -13,7 +13,7 @@ import pyproj import pyvista as pv import rasterio -from shapely.geometry import LineString, MultiLineString, MultiPolygon, Polygon, box +from shapely.geometry import LineString, MultiLineString, MultiPolygon, Point, Polygon, box from nps_active_space.utils.ais import query_ais_mxak from nps_active_space.utils.helpers import get_deployment, get_elevation, load_annotations, load_DEM, load_layered_activespace, load_studyarea from nps_active_space.scripts.run_audible_transits import AudibleTransits @@ -104,6 +104,52 @@ def active_to_linestrings(active: gpd.GeoDataFrame) -> list[LineString]: return linestrings +def vertex_z_from_coord(coord: tuple | np.ndarray) -> float: + """Return vertex elevation; missing or NaN z is treated as sea level (0 m).""" + if len(coord) < 3: + return 0.0 + z = float(coord[2]) + return 0.0 if not np.isfinite(z) else z + + +def is_surface_track(coords: np.ndarray) -> bool: + """True when every vertex is at or below sea level (vessel / surface transit).""" + return all(vertex_z_from_coord(c) <= 0.0 for c in coords) + + +def safe_dem_elevation(dem, lon: float, lat: float) -> float: + """Sample DEM elevation, returning 0 m when the point is off-raster or nodata.""" + try: + elev = float(get_elevation(dem, lon, lat)) + except (IndexError, ValueError, TypeError): + return 0.0 + return elev if np.isfinite(elev) else 0.0 + + +def _point_as_line(point: Point, span_m: float = 10.0) -> LineString: + """Convert a single-point geometry into a short renderable segment.""" + x, y = point.x, point.y + return LineString([(x, y), (x + span_m, y)]) + + +def iter_plot_linestrings(geometry) -> list[LineString]: + """Normalize annotation geometries to LineStrings PyVista can draw.""" + if geometry is None or geometry.is_empty: + return [] + if isinstance(geometry, Point): + return [_point_as_line(geometry)] + if isinstance(geometry, LineString): + if len(geometry.coords) < 2: + return [_point_as_line(Point(geometry.coords[0]))] + return [geometry] + if isinstance(geometry, MultiLineString): + lines: list[LineString] = [] + for part in geometry.geoms: + lines.extend(iter_plot_linestrings(part)) + return lines + return [] + + def densify_linestring(linestring: LineString, step_m: float) -> LineString: """Insert vertices along a line so elevation can be sampled between sparse AIS points.""" if linestring.is_empty or linestring.length <= step_m: @@ -128,7 +174,7 @@ def sea_surface_z_profile( z_vals = [] for x, y in np.array(line.coords)[:, :2]: lon, lat = to_wgs84.transform(x, y) - dem_z = float(get_elevation(dem, lon, lat)) + dem_z = safe_dem_elevation(dem, lon, lat) # NMSIM DEM uses 0 m over water; keep tracks slightly above the mesh to avoid z-fighting. z_vals.append(max(dem_z, 0.0) + offset_m) return line, np.asarray(z_vals) @@ -150,7 +196,7 @@ def create_polyline_3d( elif coords.shape[1] == 2: coords = np.column_stack((xy, np.zeros(coords.shape[0]))) assert coords.shape[1] == 3 - coords[:, 2] = np.clip(coords[:, 2], 0, 10000) # reduce magnitude of erroneous elevations + coords[:, 2] = np.clip(np.nan_to_num(coords[:, 2], nan=0.0), 0, 10000) n_points = coords.shape[0] lines = np.hstack([[n_points], np.arange(n_points)]) poly = pv.PolyData(coords) @@ -398,33 +444,27 @@ def toggle(flag): color_on=self.activespace_color ) - def _annotation_z_values(self, linestring, offset_m: float = 2.0) -> np.ndarray: - """Sample DEM elevation when track z is at or below the surface (e.g. vessel annotations).""" - coords = np.array(linestring.coords) - vertex_z = [coord[2] if len(coord) > 2 else 0.0 for coord in coords] - if all(z <= 0 for z in vertex_z): - _, z_vals = sea_surface_z_profile( - linestring, - self.dem, - self.crs, - offset_m=offset_m, - densify_step_m=self.sea_surface_densify_step_m, - ) - return z_vals - + def _annotation_z_values(self, linestring: LineString, offset_m: float = 2.0) -> np.ndarray: + """Per-vertex z for airborne tracks that may pass below the local DEM surface.""" to_wgs84 = pyproj.Transformer.from_crs(self.crs, "epsg:4326", always_xy=True) z_vals = [] - for coord in coords: + for coord in np.array(linestring.coords): x, y = coord[0], coord[1] - z = coord[2] if len(coord) > 2 else 0.0 + z = vertex_z_from_coord(coord) lon, lat = to_wgs84.transform(x, y) - dem_z = float(get_elevation(self.dem, lon, lat)) + dem_z = safe_dem_elevation(self.dem, lon, lat) if z <= 0 or z < dem_z: z_vals.append(max(dem_z, 0.0) + offset_m) else: z_vals.append(z) return np.asarray(z_vals) + def _annotation_polyline(self, linestring: LineString) -> pv.PolyData: + """Build a 3D polyline for one annotation segment.""" + if is_surface_track(np.array(linestring.coords)): + return self._sea_surface_polyline(linestring) + return create_polyline_3d(linestring, z=self._annotation_z_values(linestring)) + def _sea_surface_polyline(self, linestring: LineString) -> pv.PolyData: """Build a 3D polyline that follows the local water/ground surface.""" line, z_vals = sea_surface_z_profile( @@ -458,27 +498,36 @@ def plot_annotations(self, annotation_file: str | None = None) -> None: print(f"Sampled tracks, showing {len(selected_track_ids)} tracks = {len(annotations)} annotated segments") # clip to the area - using a box() plays more nicely with z values than the study area itself + n_loaded = len(annotations) annotations = annotations.to_crs(self.crs) - annotations = annotations.clip(box(*self.study_area.total_bounds)).explode() + annotations = annotations[~annotations.geometry.is_empty] + annotations = annotations.clip(box(*self.study_area.total_bounds)) + annotations = annotations[~annotations.geometry.is_empty].explode(ignore_index=True) # create and plot segments audible_actors = [] inaudible_actors = [] + n_plotted = 0 for _, annot in annotations.iterrows(): color = "deepskyblue" if annot["audible"] else "red" - geometry = annot["geometry"] - coords = np.array(geometry.coords) - vertex_z = [c[2] if len(c) > 2 else 0.0 for c in coords] - if all(z <= 0 for z in vertex_z): - polyline = self._sea_surface_polyline(geometry) - else: - z_vals = self._annotation_z_values(geometry) - polyline = create_polyline_3d(geometry, z=z_vals) - actor = self.plotter.add_mesh(polyline, color=color, point_size=2, line_width=2) - if annot["audible"]: - audible_actors.append(actor) - else: - inaudible_actors.append(actor) + for line in iter_plot_linestrings(annot["geometry"]): + polyline = self._annotation_polyline(line) + if polyline.n_points < 2: + continue + actor = self.plotter.add_mesh(polyline, color=color, point_size=2, line_width=2) + n_plotted += 1 + if annot["audible"]: + audible_actors.append(actor) + else: + inaudible_actors.append(actor) + + if n_plotted == 0: + print( + f"No annotation segments in view after clip " + f"({n_loaded} loaded from file; check CRS and study-area overlap)." + ) + return + print(f"Plotted {n_plotted} annotation segments") def toggle_audible(flag): for actor in audible_actors: From 01725b41ed1a5e08b41e8332a54903549f649cbc Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 11:59:22 -0800 Subject: [PATCH 05/17] more logging --- nps_active_space/scripts/viz.py | 101 ++++++++++++++++++++++++++------ tests/scripts/test_viz.py | 67 ++++++++++++++++++++- 2 files changed, 148 insertions(+), 20 deletions(-) diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index 70802b3..86346ee 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -13,9 +13,18 @@ import pyproj import pyvista as pv import rasterio +from tqdm import tqdm from shapely.geometry import LineString, MultiLineString, MultiPolygon, Point, Polygon, box from nps_active_space.utils.ais import query_ais_mxak -from nps_active_space.utils.helpers import get_deployment, get_elevation, load_annotations, load_DEM, load_layered_activespace, load_studyarea +from nps_active_space.utils.helpers import ( + get_deployment, + get_elevation, + get_logger, + load_annotations, + load_DEM, + load_layered_activespace, + load_studyarea, +) from nps_active_space.scripts.run_audible_transits import AudibleTransits from nps_active_space.utils.models import Annotations, Tracks from nps_active_space.utils.computation import NMSIM_bbox_utm @@ -150,6 +159,31 @@ def iter_plot_linestrings(geometry) -> list[LineString]: return [] +def format_annotation_summary(annotations: gpd.GeoDataFrame) -> str: + """One-line summary of a loaded annotations table for logging.""" + if annotations.empty: + return "0 segments" + n_surface = 0 + n_elevated = 0 + for geom in annotations.geometry: + if geom is None or geom.is_empty: + continue + for line in iter_plot_linestrings(geom): + if is_surface_track(np.array(line.coords)): + n_surface += 1 + else: + n_elevated += 1 + geom_types = ", ".join( + f"{k}={v}" for k, v in annotations.geometry.geom_type.value_counts().items() + ) + n_audible = int(annotations["audible"].sum()) if "audible" in annotations.columns else 0 + return ( + f"{len(annotations)} segments, {annotations['_id'].nunique()} tracks, " + f"{n_audible} audible rows, {n_surface} sea-surface / {n_elevated} elevated line(s), " + f"CRS={annotations.crs}, types: {geom_types}" + ) + + def densify_linestring(linestring: LineString, step_m: float) -> LineString: """Insert vertices along a line so elevation can be sampled between sparse AIS points.""" if linestring.is_empty or linestring.length <= step_m: @@ -256,6 +290,7 @@ def __init__( self.project_dir = cfg.read("project", "dir") self.fill_layers = fill_layers self.max_tracks = max_tracks + self.logger = get_logger("VIZ") # study area and crs self.study_area = load_studyarea(self.project_dir, self.unit, self.site, self.year) @@ -477,44 +512,67 @@ def _sea_surface_polyline(self, linestring: LineString) -> pv.PolyData: return create_polyline_3d(line, z=z_vals) def plot_annotations(self, annotation_file: str | None = None) -> None: - # load annotations if annotation_file is None: - annotations = load_annotations(self.project_dir, self.unit, self.site, self.year, only_valid=True) + self.logger.info( + f"Loading annotations from project dir for {self.deployment} (valid only)" + ) + annotations = load_annotations( + self.project_dir, self.unit, self.site, self.year, only_valid=True + ) else: - print(f"Loading annotations from custom file: {annotation_file}") + self.logger.info(f"Loading annotations from {annotation_file} (valid only)") annotations = Annotations(annotation_file, only_valid=True) - + + self.logger.info(f"Parsed annotations: {format_annotation_summary(annotations)}") if annotations.empty: - print("No annotations found, skipping.") + self.logger.info("No annotations found, skipping.") return - - # downsample n tracks if too many + track_ids = annotations["_id"].drop_duplicates() - print(f"{len(track_ids)} tracks = {len(annotations)} annotated segments") if len(track_ids) > self.max_tracks: - print(f"More than {self.max_tracks}, sampling") + self.logger.info( + f"Sampling {self.max_tracks} of {len(track_ids)} tracks " + f"({len(annotations)} segments before sample)" + ) selected_track_ids = track_ids.sample(self.max_tracks, replace=False, random_state=2) annotations = annotations[annotations["_id"].isin(selected_track_ids)] - print(f"Sampled tracks, showing {len(selected_track_ids)} tracks = {len(annotations)} annotated segments") + self.logger.info( + f"After sample: {format_annotation_summary(annotations)}" + ) - # clip to the area - using a box() plays more nicely with z values than the study area itself n_loaded = len(annotations) + self.logger.info(f"Reprojecting annotations to {self.crs} and clipping to study area") annotations = annotations.to_crs(self.crs) + n_empty = int(annotations.geometry.is_empty.sum()) + if n_empty: + self.logger.info(f"Dropping {n_empty} empty geometries before clip") annotations = annotations[~annotations.geometry.is_empty] annotations = annotations.clip(box(*self.study_area.total_bounds)) annotations = annotations[~annotations.geometry.is_empty].explode(ignore_index=True) + self.logger.info( + f"After clip: {len(annotations)} segments " + f"({n_loaded - len(annotations)} removed outside study area)" + ) - # create and plot segments audible_actors = [] inaudible_actors = [] n_plotted = 0 - for _, annot in annotations.iterrows(): + n_skipped = 0 + for _, annot in tqdm( + annotations.iterrows(), + total=len(annotations), + desc="Plotting annotations", + unit="segment", + ): color = "deepskyblue" if annot["audible"] else "red" for line in iter_plot_linestrings(annot["geometry"]): polyline = self._annotation_polyline(line) if polyline.n_points < 2: + n_skipped += 1 continue - actor = self.plotter.add_mesh(polyline, color=color, point_size=2, line_width=2) + actor = self.plotter.add_mesh( + polyline, color=color, point_size=2, line_width=2 + ) n_plotted += 1 if annot["audible"]: audible_actors.append(actor) @@ -522,12 +580,17 @@ def plot_annotations(self, annotation_file: str | None = None) -> None: inaudible_actors.append(actor) if n_plotted == 0: - print( - f"No annotation segments in view after clip " - f"({n_loaded} loaded from file; check CRS and study-area overlap)." + self.logger.info( + f"No annotation segments plotted ({n_loaded} loaded; " + f"{n_skipped} degenerate lines skipped). " + "Check CRS and study-area overlap." ) return - print(f"Plotted {n_plotted} annotation segments") + self.logger.info( + f"Plotted {n_plotted} annotation line(s) " + f"({len(audible_actors)} audible, {len(inaudible_actors)} inaudible; " + f"{n_skipped} skipped)" + ) def toggle_audible(flag): for actor in audible_actors: diff --git a/tests/scripts/test_viz.py b/tests/scripts/test_viz.py index 1b2d637..16ef331 100644 --- a/tests/scripts/test_viz.py +++ b/tests/scripts/test_viz.py @@ -3,17 +3,21 @@ import numpy as np import pytest -from shapely.geometry import LineString, Point +from shapely.geometry import LineString, MultiLineString, Point from nps_active_space.scripts.viz import ( create_polyline_3d, densify_linestring, + format_annotation_summary, + is_surface_track, + iter_plot_linestrings, parse_deployment, parse_existing_file, parse_iso_date, parse_max_tracks, sea_surface_z_profile, track_points_to_linestring, + vertex_z_from_coord, ) @@ -123,7 +127,68 @@ def test_negative_dem_clamped_to_sea_level(self, monkeypatch): assert z_vals.tolist() == [2.0, 2.0] +class TestFormatAnnotationSummary: + def test_empty(self): + import geopandas as gpd + + assert format_annotation_summary(gpd.GeoDataFrame()) == "0 segments" + + def test_includes_track_and_surface_counts(self): + import geopandas as gpd + + gdf = gpd.GeoDataFrame( + { + "_id": ["T1", "T1"], + "audible": [True, False], + }, + geometry=[ + LineString([(-136.0, 58.0), (-136.1, 58.1)]), + LineString([(0.0, 0.0, 100.0), (1.0, 1.0, 200.0)]), + ], + crs="EPSG:4326", + ) + summary = format_annotation_summary(gdf) + assert "2 segments" in summary + assert "1 tracks" in summary + assert "1 audible" in summary + assert "1 sea-surface" in summary + assert "1 elevated" in summary + + +class TestVertexZAndSurfaceTrack: + def test_missing_z_is_sea_level(self): + assert vertex_z_from_coord((-136.0, 58.0)) == 0.0 + + def test_nan_z_is_sea_level(self): + assert vertex_z_from_coord((-136.0, 58.0, float("nan"))) == 0.0 + + def test_nan_z_counts_as_surface_track(self): + line = LineString([(-136.0, 58.0, float("nan")), (-136.1, 58.1, float("nan"))]) + assert is_surface_track(np.array(line.coords)) + + def test_positive_z_is_not_surface_track(self): + line = LineString([(0, 0, 100.0), (1, 1, 200.0)]) + assert not is_surface_track(np.array(line.coords)) + + +class TestIterPlotLinestrings: + def test_point_becomes_short_line(self): + lines = iter_plot_linestrings(Point(1.0, 2.0)) + assert len(lines) == 1 + assert len(lines[0].coords) == 2 + + def test_multilinestring_splits(self): + geom = MultiLineString([LineString([(0, 0), (1, 1)]), LineString([(2, 2), (3, 3)])]) + assert len(iter_plot_linestrings(geom)) == 2 + + class TestCreatePolyline3d: + def test_nan_z_becomes_zero(self): + line = LineString([(0, 0, float("nan")), (1, 0, float("nan"))]) + poly = create_polyline_3d(line, z=np.array([float("nan"), 5.0])) + assert poly.points[0, 2] == pytest.approx(0.0) + assert poly.points[1, 2] == pytest.approx(5.0) + def test_rejects_mismatched_z_length(self): line = LineString([(0, 0), (1, 1), (2, 2)]) with pytest.raises(ValueError, match="one value per vertex"): From 718d045ec9c6868f2396b25bea72f41b07162719 Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 12:13:54 -0800 Subject: [PATCH 06/17] more twaeks + script for downsmapling geojson --- nps_active_space/ground_truthing/segments.py | 72 ++++++++++++- .../scripts/downsample_annotations.py | 85 +++++++++++++++ nps_active_space/scripts/viz.py | 100 ++++++++++++------ tests/ground_truthing/test_segments.py | 36 +++++++ tests/scripts/test_viz.py | 15 +++ 5 files changed, 271 insertions(+), 37 deletions(-) create mode 100644 nps_active_space/scripts/downsample_annotations.py diff --git a/nps_active_space/ground_truthing/segments.py b/nps_active_space/ground_truthing/segments.py index 3c36a08..0a81a5f 100644 --- a/nps_active_space/ground_truthing/segments.py +++ b/nps_active_space/ground_truthing/segments.py @@ -1,11 +1,65 @@ import datetime as dt import geopandas as gpd +import numpy as np from shapely.geometry import LineString, Point # Mutable [start, end] pair on the time_audible axis (lists are updated in place). AudibleRange = list[dt.datetime] +# Splines are 1 Hz; storing every vertex makes GeoJSON huge. Reload uses start_dt/end_dt only. +ANNOTATION_MAX_VERTICES = 64 +COORD_DECIMALS = 6 +Z_DECIMALS = 1 + + +def _round_coord(coord: tuple) -> tuple: + """Round WGS84 xy and elevation for compact GeoJSON.""" + if len(coord) < 3: + return (round(coord[0], COORD_DECIMALS), round(coord[1], COORD_DECIMALS)) + return ( + round(coord[0], COORD_DECIMALS), + round(coord[1], COORD_DECIMALS), + round(coord[2], Z_DECIMALS), + ) + + +def downsample_linestring( + line: LineString, + max_vertices: int = ANNOTATION_MAX_VERTICES, +) -> LineString: + """Reduce vertex count for annotation storage and viz (keeps endpoints).""" + if line.is_empty: + return line + coords = list(line.coords) + if len(coords) <= max_vertices: + return LineString([_round_coord(c) for c in coords]) + indices = np.unique(np.linspace(0, len(coords) - 1, max_vertices, dtype=int)) + sampled = [_round_coord(coords[i]) for i in indices] + return LineString(sampled) + + +def storage_geometry_from_points(geometries: list) -> Point | LineString: + """Build a compact Point or downsampled LineString for GeoJSON export.""" + if len(geometries) == 1: + geom = geometries[0] + if isinstance(geom, Point): + return Point(_round_coord(geom.coords[0])) + return downsample_linestring(LineString([geom])) + return downsample_linestring(LineString(geometries)) + + +def compact_geometry( + geom: Point | LineString, + max_vertices: int = ANNOTATION_MAX_VERTICES, +) -> Point | LineString: + """Downsample and round a stored annotation geometry.""" + if isinstance(geom, Point): + return Point(_round_coord(geom.coords[0])) + if isinstance(geom, LineString): + return downsample_linestring(geom, max_vertices=max_vertices) + return geom + def collapse_audible_ranges(ranges: list[AudibleRange]) -> list[AudibleRange]: """Collapse overlapping audible ranges into single audible ranges. @@ -85,8 +139,10 @@ def build_annotation_segments( 'valid': valid, 'audible': False, 'note': note, - 'geometry': points.geometry.iat[0] if points.shape[0] == 1 - else LineString(points.geometry.tolist()) + 'geometry': storage_geometry_from_points( + [points.geometry.iat[0]] if points.shape[0] == 1 + else list(points.geometry) + ) }) else: @@ -112,7 +168,9 @@ def build_annotation_segments( 'valid': True, 'audible': False, 'note': note, - 'geometry': LineString(first_inaudible_segment.geometry.tolist())} + 'geometry': storage_geometry_from_points( + list(first_inaudible_segment.geometry) + )} ) # add segments for each audible range and the inaudible range following it @@ -126,7 +184,9 @@ def build_annotation_segments( 'valid': True, 'audible': True, 'note': note, - 'geometry': LineString(audible_segment.geometry.tolist())} + 'geometry': storage_geometry_from_points( + list(audible_segment.geometry) + )} ) if i+1 < len(audible_ranges): @@ -143,7 +203,9 @@ def build_annotation_segments( 'valid': True, 'audible': False, 'note': note, - 'geometry': LineString(inaudible_segment.geometry.tolist())} + 'geometry': storage_geometry_from_points( + list(inaudible_segment.geometry) + )} ) gdf = gpd.GeoDataFrame(segments, geometry='geometry', crs=points.crs) diff --git a/nps_active_space/scripts/downsample_annotations.py b/nps_active_space/scripts/downsample_annotations.py new file mode 100644 index 0000000..337ac5d --- /dev/null +++ b/nps_active_space/scripts/downsample_annotations.py @@ -0,0 +1,85 @@ +"""Downsample dense annotation GeoJSON (1 Hz splines) for smaller files.""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +import geopandas as gpd +from shapely.geometry import LineString + +from nps_active_space.ground_truthing.segments import ( + ANNOTATION_MAX_VERTICES, + compact_geometry, +) + + +def _vertex_count(geom) -> int: + if isinstance(geom, LineString): + return len(geom.coords) + return 1 + + +def downsample_annotations_file( + input_path: Path, + output_path: Path, + *, + max_vertices: int = ANNOTATION_MAX_VERTICES, +) -> tuple[int, int]: + """Read annotation GeoJSON, compact geometries, write output. Returns (before, after) vertex totals.""" + gdf = gpd.read_file(input_path) + before = sum(_vertex_count(g) for g in gdf.geometry) + gdf = gdf.copy() + gdf["geometry"] = gdf.geometry.map( + lambda g: compact_geometry(g, max_vertices=max_vertices) + ) + after = sum(_vertex_count(g) for g in gdf.geometry) + output_path.parent.mkdir(parents=True, exist_ok=True) + gdf.to_file(output_path, driver="GeoJSON") + return before, after + + +def main() -> None: + parser = argparse.ArgumentParser( + description="Downsample dense ground-truthing annotation GeoJSON." + ) + parser.add_argument("input", type=Path, help="Input .geojson") + parser.add_argument( + "output", + type=Path, + nargs="?", + help="Output .geojson (default: _compact.geojson)", + ) + parser.add_argument( + "--max-vertices", + type=int, + default=ANNOTATION_MAX_VERTICES, + help=f"Max vertices per segment (default: {ANNOTATION_MAX_VERTICES})", + ) + args = parser.parse_args() + + input_path = args.input.resolve() + if not input_path.is_file(): + parser.error(f"input not found: {input_path}") + + output_path = ( + args.output.resolve() + if args.output is not None + else input_path.with_name(f"{input_path.stem}_compact.geojson") + ) + + before_bytes = input_path.stat().st_size + before_verts, after_verts = downsample_annotations_file( + input_path, + output_path, + max_vertices=args.max_vertices, + ) + after_bytes = output_path.stat().st_size + + print(f"Wrote {output_path}") + print(f" vertices: {before_verts:,} -> {after_verts:,}") + print(f" file size: {before_bytes / 1e6:.1f} MB -> {after_bytes / 1e6:.1f} MB") + + +if __name__ == "__main__": + main() diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index 86346ee..d448301 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -85,6 +85,24 @@ def parse_max_tracks(value: str) -> int: return n +def resolve_viz_plot_flags( + *, + active_space: bool = False, + annotations: bool = False, + audible_transits: bool = False, + vessels: bool = False, + plot_all: bool = False, + annotation_file: str | None = None, + transits_pkl: str | None = None, +) -> tuple[bool, bool, bool, bool]: + """Map CLI flags to Visualizer layer toggles.""" + do_active = active_space or plot_all + do_annotations = annotations or plot_all or annotation_file is not None + do_transits = audible_transits or plot_all or transits_pkl is not None + do_vessels = vessels + return do_active, do_annotations, do_transits, do_vessels + + # helper functions ============================================ def polygon_to_mesh(polygon: Polygon, z: float) -> pv.PolyData: @@ -290,7 +308,7 @@ def __init__( self.project_dir = cfg.read("project", "dir") self.fill_layers = fill_layers self.max_tracks = max_tracks - self.logger = get_logger("VIZ") + self.logger = get_logger("VIZ", verbose=True) # study area and crs self.study_area = load_studyarea(self.project_dir, self.unit, self.site, self.year) @@ -311,13 +329,30 @@ def __init__( self.plot_vessel_tracks(vessel_start_date, vessel_end_date) # configure plot parameters and display - self.plotter.camera_position = "xz" - self.plotter.camera.elevation = 45 self.plotter.enable_terrain_style() self.setup_z_scale() self.plotter.add_title(f"{unit}{site}{year}", font_size=12) self.plotter.show_axes() + self.plotter.reset_camera() + self.plotter.camera.elevation = 30 self.plotter.show() + + def _status(self, message: str) -> None: + """Print and log a user-facing status line (visible when run via python).""" + print(message, flush=True) + self.logger.info(message) + + def _add_track_line( + self, polyline: pv.PolyData, *, color: str, line_width: int = 4 + ): + """Add a track/annotation polyline that reads clearly over the DEM.""" + return self.plotter.add_mesh( + polyline, + color=color, + line_width=line_width, + render_lines_as_tubes=True, + point_size=4, + ) def plot_dem(self, show_scalar_bar: bool = False) -> None: # load DEM @@ -513,43 +548,46 @@ def _sea_surface_polyline(self, linestring: LineString) -> pv.PolyData: def plot_annotations(self, annotation_file: str | None = None) -> None: if annotation_file is None: - self.logger.info( + self._status( f"Loading annotations from project dir for {self.deployment} (valid only)" ) annotations = load_annotations( self.project_dir, self.unit, self.site, self.year, only_valid=True ) else: - self.logger.info(f"Loading annotations from {annotation_file} (valid only)") + self._status(f"Loading annotations from {annotation_file} (valid only)") annotations = Annotations(annotation_file, only_valid=True) - self.logger.info(f"Parsed annotations: {format_annotation_summary(annotations)}") + self._status(f"Parsed annotations: {format_annotation_summary(annotations)}") if annotations.empty: - self.logger.info("No annotations found, skipping.") + self._status("No annotations found, skipping.") return track_ids = annotations["_id"].drop_duplicates() if len(track_ids) > self.max_tracks: - self.logger.info( + self._status( f"Sampling {self.max_tracks} of {len(track_ids)} tracks " f"({len(annotations)} segments before sample)" ) selected_track_ids = track_ids.sample(self.max_tracks, replace=False, random_state=2) annotations = annotations[annotations["_id"].isin(selected_track_ids)] - self.logger.info( - f"After sample: {format_annotation_summary(annotations)}" - ) + self._status(f"After sample: {format_annotation_summary(annotations)}") n_loaded = len(annotations) - self.logger.info(f"Reprojecting annotations to {self.crs} and clipping to study area") + self._status(f"Reprojecting annotations to {self.crs} and clipping to study area") annotations = annotations.to_crs(self.crs) + ann_bounds = annotations.total_bounds + study_bounds = self.study_area.total_bounds + self._status(f"Annotation bounds: {ann_bounds}") + self._status(f"Study area bounds: {study_bounds}") + n_empty = int(annotations.geometry.is_empty.sum()) if n_empty: - self.logger.info(f"Dropping {n_empty} empty geometries before clip") + self._status(f"Dropping {n_empty} empty geometries before clip") annotations = annotations[~annotations.geometry.is_empty] - annotations = annotations.clip(box(*self.study_area.total_bounds)) + annotations = annotations.clip(box(*study_bounds)) annotations = annotations[~annotations.geometry.is_empty].explode(ignore_index=True) - self.logger.info( + self._status( f"After clip: {len(annotations)} segments " f"({n_loaded - len(annotations)} removed outside study area)" ) @@ -570,9 +608,7 @@ def plot_annotations(self, annotation_file: str | None = None) -> None: if polyline.n_points < 2: n_skipped += 1 continue - actor = self.plotter.add_mesh( - polyline, color=color, point_size=2, line_width=2 - ) + actor = self._add_track_line(polyline, color=color) n_plotted += 1 if annot["audible"]: audible_actors.append(actor) @@ -580,13 +616,13 @@ def plot_annotations(self, annotation_file: str | None = None) -> None: inaudible_actors.append(actor) if n_plotted == 0: - self.logger.info( + self._status( f"No annotation segments plotted ({n_loaded} loaded; " f"{n_skipped} degenerate lines skipped). " "Check CRS and study-area overlap." ) return - self.logger.info( + self._status( f"Plotted {n_plotted} annotation line(s) " f"({len(audible_actors)} audible, {len(inaudible_actors)} inaudible; " f"{n_skipped} skipped)" @@ -704,12 +740,7 @@ def plot_vessel_tracks( if line.is_empty or line.length == 0: continue polyline = self._sea_surface_polyline(line) - actor = self.plotter.add_mesh( - polyline, - color=self.vessel_track_color, - point_size=2, - line_width=3, - ) + actor = self._add_track_line(polyline, color=self.vessel_track_color, line_width=5) actors.append(actor) if not actors: @@ -783,12 +814,12 @@ def toggle_z_scale(flag): parser.add_argument( "--annotation-file", type=lambda p: parse_existing_file(p, arg_name="--annotation-file"), - help="Path to .geojson file to load annotations from, if not the default.", + help="Path to .geojson annotations (implies -a).", ) parser.add_argument( "--transits-pkl", type=lambda p: parse_existing_file(p, arg_name="--transits-pkl"), - help="Path to .pkl file to load audible transits from, if not the default.", + help="Path to .pkl audible transits (implies -t).", ) parser.add_argument( "--start-date", @@ -810,10 +841,15 @@ def toggle_z_scale(flag): unit, site, year = args.deployment print(unit, site, year) - do_active = args.active_space or args.all - do_annotations = args.annotations or args.all - do_transits = args.audible_transits or args.all - do_vessels = args.vessels + do_active, do_annotations, do_transits, do_vessels = resolve_viz_plot_flags( + active_space=args.active_space, + annotations=args.annotations, + audible_transits=args.audible_transits, + vessels=args.vessels, + plot_all=args.all, + annotation_file=args.annotation_file, + transits_pkl=args.transits_pkl, + ) Visualizer(unit, site, year, args.environment, do_active, args.gain, diff --git a/tests/ground_truthing/test_segments.py b/tests/ground_truthing/test_segments.py index a6f2448..3dfa937 100644 --- a/tests/ground_truthing/test_segments.py +++ b/tests/ground_truthing/test_segments.py @@ -1,12 +1,15 @@ import datetime as dt import pandas as pd +import pytest from pandas.testing import assert_frame_equal from shapely.geometry import LineString from nps_active_space.ground_truthing.segments import ( + ANNOTATION_MAX_VERTICES, build_annotation_segments, collapse_audible_ranges, + downsample_linestring, ) from helpers import SEGMENT_TABULAR_COLS, _tabular, make_track_points @@ -168,3 +171,36 @@ def test_note_propagates_to_segments(self): } ) assert_frame_equal(actual, expected, check_dtype=False) + + def test_long_track_geometry_downsampled(self): + points = make_track_points(500, freq="s") + t = points.point_dt + result = build_annotation_segments( + "T1", + points, + audible_ranges=[[t.iat[0], t.iat[-1]]], + ) + geom = result.iloc[0].geometry + assert isinstance(geom, LineString) + assert len(geom.coords) <= ANNOTATION_MAX_VERTICES + assert geom.coords[0][0] == pytest.approx(0.0, abs=1e-6) + assert geom.coords[0][1] == pytest.approx(0.0, abs=1e-6) + assert geom.coords[-1][0] == pytest.approx(498.0, abs=1e-6) + assert geom.coords[-1][1] == pytest.approx(498.0, abs=1e-6) + + +class TestDownsampleLinestring: + def test_caps_vertices_and_preserves_endpoints(self): + coords = [(float(i), float(i), float(i)) for i in range(1000)] + line = LineString(coords) + out = downsample_linestring(line, max_vertices=32) + assert len(out.coords) == 32 + assert out.coords[0] == pytest.approx(coords[0], abs=1e-6) + assert out.coords[-1] == pytest.approx(coords[-1], abs=1e-6) + + def test_short_line_only_rounds(self): + line = LineString([(1.123456789, 2.987654321, 100.55), (3.0, 4.0, 200.44)]) + out = downsample_linestring(line) + assert len(out.coords) == 2 + assert out.coords[0] == (1.123457, 2.987654, 100.5) + assert out.coords[1] == (3.0, 4.0, 200.4) diff --git a/tests/scripts/test_viz.py b/tests/scripts/test_viz.py index 16ef331..b4af2b6 100644 --- a/tests/scripts/test_viz.py +++ b/tests/scripts/test_viz.py @@ -15,6 +15,7 @@ parse_existing_file, parse_iso_date, parse_max_tracks, + resolve_viz_plot_flags, sea_surface_z_profile, track_points_to_linestring, vertex_z_from_coord, @@ -64,6 +65,20 @@ def test_missing_file(self, tmp_path: Path): parse_existing_file(str(tmp_path / "missing.geojson"), arg_name="--transits-pkl") +class TestResolveVizPlotFlags: + def test_annotation_file_implies_annotations(self): + flags = resolve_viz_plot_flags(annotation_file="/tmp/a.geojson") + assert flags == (False, True, False, False) + + def test_transits_pkl_implies_transits(self): + flags = resolve_viz_plot_flags(transits_pkl="/tmp/t.pkl") + assert flags == (False, False, True, False) + + def test_all_enables_standard_layers(self): + flags = resolve_viz_plot_flags(plot_all=True) + assert flags == (True, True, True, False) + + class TestParseMaxTracks: def test_valid_positive_int(self): assert parse_max_tracks("500") == 500 From 339d0326f348e3bc7ede3358659ea0ffab24602a Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 12:21:52 -0800 Subject: [PATCH 07/17] cursor perf improvmenets on dem loading --- nps_active_space/scripts/viz.py | 57 +++++++++++++++++++++++++------ nps_active_space/utils/helpers.py | 6 ++-- tests/scripts/test_viz.py | 9 +++++ 3 files changed, 58 insertions(+), 14 deletions(-) diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index d448301..68a2507 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -314,6 +314,7 @@ def __init__( self.study_area = load_studyarea(self.project_dir, self.unit, self.site, self.year) self.crs = NMSIM_bbox_utm(self.study_area) self.study_area = self.study_area.to_crs(self.crs) + self._to_wgs84 = pyproj.Transformer.from_crs(self.crs, "epsg:4326", always_xy=True) # plot each element self.plotter = pv.Plotter() @@ -353,6 +354,29 @@ def _add_track_line( render_lines_as_tubes=True, point_size=4, ) + + def _add_annotation_lines( + self, polylines: list[pv.PolyData], *, color: str, line_width: int = 3 + ): + """Add many annotation segments as one mesh (much faster than per-segment tubes).""" + polylines = [p for p in polylines if p.n_points >= 2] + if not polylines: + return None + mesh = polylines[0] if len(polylines) == 1 else pv.merge(polylines) + return self.plotter.add_mesh( + mesh, + color=color, + line_width=line_width, + render_lines_as_tubes=False, + ) + + @staticmethod + def _flat_sea_surface_polyline( + linestring: LineString, offset_m: float + ) -> pv.PolyData: + """Sea-level track at a constant offset — no DEM resampling.""" + z = np.full(len(linestring.coords), offset_m) + return create_polyline_3d(linestring, z=z) def plot_dem(self, show_scalar_bar: bool = False) -> None: # load DEM @@ -516,12 +540,11 @@ def toggle(flag): def _annotation_z_values(self, linestring: LineString, offset_m: float = 2.0) -> np.ndarray: """Per-vertex z for airborne tracks that may pass below the local DEM surface.""" - to_wgs84 = pyproj.Transformer.from_crs(self.crs, "epsg:4326", always_xy=True) z_vals = [] for coord in np.array(linestring.coords): x, y = coord[0], coord[1] z = vertex_z_from_coord(coord) - lon, lat = to_wgs84.transform(x, y) + lon, lat = self._to_wgs84.transform(x, y) dem_z = safe_dem_elevation(self.dem, lon, lat) if z <= 0 or z < dem_z: z_vals.append(max(dem_z, 0.0) + offset_m) @@ -532,11 +555,13 @@ def _annotation_z_values(self, linestring: LineString, offset_m: float = 2.0) -> def _annotation_polyline(self, linestring: LineString) -> pv.PolyData: """Build a 3D polyline for one annotation segment.""" if is_surface_track(np.array(linestring.coords)): - return self._sea_surface_polyline(linestring) + return self._flat_sea_surface_polyline(linestring, self.sea_surface_offset_m) return create_polyline_3d(linestring, z=self._annotation_z_values(linestring)) def _sea_surface_polyline(self, linestring: LineString) -> pv.PolyData: """Build a 3D polyline that follows the local water/ground surface.""" + if is_surface_track(np.array(linestring.coords)): + return self._flat_sea_surface_polyline(linestring, self.sea_surface_offset_m) line, z_vals = sea_surface_z_profile( linestring, self.dem, @@ -594,26 +619,37 @@ def plot_annotations(self, annotation_file: str | None = None) -> None: audible_actors = [] inaudible_actors = [] + audible_polylines: list[pv.PolyData] = [] + inaudible_polylines: list[pv.PolyData] = [] n_plotted = 0 n_skipped = 0 for _, annot in tqdm( annotations.iterrows(), total=len(annotations), - desc="Plotting annotations", + desc="Building annotation lines", unit="segment", ): - color = "deepskyblue" if annot["audible"] else "red" for line in iter_plot_linestrings(annot["geometry"]): polyline = self._annotation_polyline(line) if polyline.n_points < 2: n_skipped += 1 continue - actor = self._add_track_line(polyline, color=color) n_plotted += 1 if annot["audible"]: - audible_actors.append(actor) + audible_polylines.append(polyline) else: - inaudible_actors.append(actor) + inaudible_polylines.append(polyline) + + audible_actor = self._add_annotation_lines( + audible_polylines, color=self.audible_annotation_color + ) + if audible_actor is not None: + audible_actors.append(audible_actor) + inaudible_actor = self._add_annotation_lines( + inaudible_polylines, color=self.inaudible_annotation_color + ) + if inaudible_actor is not None: + inaudible_actors.append(inaudible_actor) if n_plotted == 0: self._status( @@ -623,8 +659,9 @@ def plot_annotations(self, annotation_file: str | None = None) -> None: ) return self._status( - f"Plotted {n_plotted} annotation line(s) " - f"({len(audible_actors)} audible, {len(inaudible_actors)} inaudible; " + f"Plotted {n_plotted} annotation line(s) in " + f"{len(audible_actors) + len(inaudible_actors)} mesh(es) " + f"({len(audible_polylines)} audible, {len(inaudible_polylines)} inaudible; " f"{n_skipped} skipped)" ) diff --git a/nps_active_space/utils/helpers.py b/nps_active_space/utils/helpers.py index 1599650..f46a8bd 100644 --- a/nps_active_space/utils/helpers.py +++ b/nps_active_space/utils/helpers.py @@ -161,13 +161,11 @@ def load_DEM(project_dir: str, unit: str, site: str): def get_elevation(DEM, lon: float, lat: float) -> float: - """Read elevation at a certain lat/lon from a DEM. - """ + """Read elevation at a certain lat/lon from a DEM.""" proj = Transformer.from_crs("epsg:4326", DEM.crs, always_xy=True) x, y = proj.transform(lon, lat) row, col = DEM.index(x, y) - elevation = DEM.read(1)[row, col] - return elevation + return float(DEM.read(1, window=((row, row + 1), (col, col + 1)))[0, 0]) def load_studyarea(project_dir: str, unit: str, site: str, year: int, crs: str = None): diff --git a/tests/scripts/test_viz.py b/tests/scripts/test_viz.py index b4af2b6..dc96d64 100644 --- a/tests/scripts/test_viz.py +++ b/tests/scripts/test_viz.py @@ -6,6 +6,7 @@ from shapely.geometry import LineString, MultiLineString, Point from nps_active_space.scripts.viz import ( + Visualizer, create_polyline_3d, densify_linestring, format_annotation_summary, @@ -170,6 +171,14 @@ def test_includes_track_and_surface_counts(self): assert "1 elevated" in summary +class TestFlatSeaSurfacePolyline: + def test_uses_constant_offset_without_dem(self): + line = LineString([(0, 0), (1000, 0), (2000, 0)]) + poly = Visualizer._flat_sea_surface_polyline(line, offset_m=5.0) + assert poly.n_points == 3 + assert np.all(poly.points[:, 2] == pytest.approx(5.0)) + + class TestVertexZAndSurfaceTrack: def test_missing_z_is_sea_level(self): assert vertex_z_from_coord((-136.0, 58.0)) == 0.0 From 4cf89e9e02a400885de3ebc4fa594cf41d4890b4 Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 12:26:45 -0800 Subject: [PATCH 08/17] compass --- nps_active_space/scripts/viz.py | 18 +++++++++++++++++- tests/scripts/test_viz.py | 10 ++++++++++ 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index 68a2507..6544387 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -279,6 +279,8 @@ class Visualizer(): z_scale_toggle_color = "black" # button toggling z scale sea_surface_offset_m = 5.0 sea_surface_densify_step_m = 100.0 + axes_viewport = (0.0, 0.0, 0.14, 0.14) + compass_viewport = (0.14, 0.0, 0.26, 0.14) def __init__( self, @@ -333,7 +335,7 @@ def __init__( self.plotter.enable_terrain_style() self.setup_z_scale() self.plotter.add_title(f"{unit}{site}{year}", font_size=12) - self.plotter.show_axes() + self.setup_orientation_widgets() self.plotter.reset_camera() self.plotter.camera.elevation = 30 self.plotter.show() @@ -796,6 +798,20 @@ def toggle(flag): color_on=self.vessel_track_color, ) + def setup_orientation_widgets(self) -> None: + """Bottom-left XYZ axes and geographic north compass (UTM +Y).""" + self.plotter.add_axes( + viewport=self.axes_viewport, + interactive=False, + line_width=2, + ) + self.plotter.add_north_arrow_widget( + viewport=self.compass_viewport, + interactive=False, + color="#2f4f4f", + line_width=2, + ) + def setup_z_scale(self) -> None: self.plotter.set_scale(1, 1, 2) # easier to make out elevation differences diff --git a/tests/scripts/test_viz.py b/tests/scripts/test_viz.py index dc96d64..aa20580 100644 --- a/tests/scripts/test_viz.py +++ b/tests/scripts/test_viz.py @@ -171,6 +171,16 @@ def test_includes_track_and_surface_counts(self): assert "1 elevated" in summary +class TestOrientationWidgets: + def test_compass_viewport_sits_right_of_axes(self): + ax = Visualizer.axes_viewport + compass = Visualizer.compass_viewport + assert ax[0] == 0.0 + assert compass[0] == pytest.approx(ax[2]) + assert compass[2] > compass[0] + assert compass[3] == pytest.approx(ax[3]) + + class TestFlatSeaSurfacePolyline: def test_uses_constant_offset_without_dem(self): line = LineString([(0, 0), (1000, 0), (2000, 0)]) From ccaf6e796ae48fbcb36b0179a2c1e7836d06aa2e Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 12:34:11 -0800 Subject: [PATCH 09/17] test --- nps_active_space/scripts/viz.py | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index 6544387..d926019 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -103,6 +103,22 @@ def resolve_viz_plot_flags( return do_active, do_annotations, do_transits, do_vessels +def utm_orientation_axes_kwargs() -> dict: + """Orientation-marker kwargs for UTM scenes (+X east, +Y north). + + PyVista/VTK allow only one vtkOrientationMarkerWidget per renderer, so axes + and north-arrow widgets cannot coexist as separate corner widgets. + """ + return { + "interactive": False, + "line_width": 2, + "xlabel": "E", + "ylabel": "N", + "zlabel": "Z", + "viewport": (0, 0, 0.2, 0.2), + } + + # helper functions ============================================ def polygon_to_mesh(polygon: Polygon, z: float) -> pv.PolyData: @@ -279,8 +295,6 @@ class Visualizer(): z_scale_toggle_color = "black" # button toggling z scale sea_surface_offset_m = 5.0 sea_surface_densify_step_m = 100.0 - axes_viewport = (0.0, 0.0, 0.14, 0.14) - compass_viewport = (0.14, 0.0, 0.26, 0.14) def __init__( self, @@ -799,18 +813,8 @@ def toggle(flag): ) def setup_orientation_widgets(self) -> None: - """Bottom-left XYZ axes and geographic north compass (UTM +Y).""" - self.plotter.add_axes( - viewport=self.axes_viewport, - interactive=False, - line_width=2, - ) - self.plotter.add_north_arrow_widget( - viewport=self.compass_viewport, - interactive=False, - color="#2f4f4f", - line_width=2, - ) + """Bottom-left E/N/Z axes (+Y = north in UTM).""" + self.plotter.add_axes(**utm_orientation_axes_kwargs()) def setup_z_scale(self) -> None: self.plotter.set_scale(1, 1, 2) # easier to make out elevation differences From e08442ef0fc3b7dca31ff5c5532d47b9fe33b20f Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 12:37:58 -0800 Subject: [PATCH 10/17] improve aircraft loading --- nps_active_space/scripts/viz.py | 94 +++++++++++++++++++++++++++------ tests/scripts/test_viz.py | 38 ++++++++++--- 2 files changed, 109 insertions(+), 23 deletions(-) diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index d926019..81e5379 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -160,6 +160,77 @@ def is_surface_track(coords: np.ndarray) -> bool: return all(vertex_z_from_coord(c) <= 0.0 for c in coords) +def is_airborne_track(coords: np.ndarray) -> bool: + """True when every vertex has positive stored elevation (aircraft spline z).""" + return all(vertex_z_from_coord(c) > 0.0 for c in coords) + + +def stored_vertex_z(linestring: LineString) -> np.ndarray: + """Elevation from annotation geometry coordinates.""" + return np.array([vertex_z_from_coord(c) for c in linestring.coords]) + + +class DemElevationSampler: + """Sample elevations from an in-memory DEM band (no per-point raster I/O).""" + + def __init__(self, dem, band: np.ndarray, plot_crs: str) -> None: + self.dem = dem + self.band = np.asarray(band, dtype=float).copy() + if dem.nodata is not None: + self.band[self.band == dem.nodata] = 0.0 + self.band[self.band > 9000] = 0.0 + self._to_dem = pyproj.Transformer.from_crs(plot_crs, dem.crs, always_xy=True) + + def sample_utm_many(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + dem_x, dem_y = self._to_dem.transform(x, y) + rows, cols = self.dem.index(dem_x, dem_y) + rows = np.atleast_1d(rows).astype(int) + cols = np.atleast_1d(cols).astype(int) + in_bounds = ( + (rows >= 0) + & (rows < self.band.shape[0]) + & (cols >= 0) + & (cols < self.band.shape[1]) + ) + elev = np.zeros(rows.shape[0], dtype=float) + if in_bounds.any(): + samples = self.band[rows[in_bounds], cols[in_bounds]] + elev[in_bounds] = np.where(np.isfinite(samples), samples, 0.0) + return elev + + +def annotation_z_profile( + linestring: LineString, + dem_sampler: DemElevationSampler | None, + *, + offset_m: float = 2.0, +) -> np.ndarray: + """Per-vertex z for annotation segments. + + Aircraft splines already carry MSL elevation in geometry; only vertices at or + below sea level need a local DEM clamp (vessels use the flat sea-surface path). + """ + z_stored = stored_vertex_z(linestring) + if np.all(z_stored > 0): + return z_stored + if dem_sampler is None: + return z_stored + coords = np.array(linestring.coords) + need_dem = z_stored <= 0 + if not need_dem.any(): + return z_stored + dem_z = dem_sampler.sample_utm_many(coords[need_dem, 0], coords[need_dem, 1]) + z_out = z_stored.copy() + for j, i in enumerate(np.where(need_dem)[0]): + z = z_stored[i] + dz = float(dem_z[j]) + if z <= 0 or z < dz: + z_out[i] = max(dz, 0.0) + offset_m + return z_out + + def safe_dem_elevation(dem, lon: float, lat: float) -> float: """Sample DEM elevation, returning 0 m when the point is off-raster or nodata.""" try: @@ -402,6 +473,7 @@ def plot_dem(self, show_scalar_bar: bool = False) -> None: if dem.nodata is not None: data[data == dem.nodata] = 0 data[data > 9000] = 0 # higher than any elevation on earth, should be nodata + self._dem_sampler = DemElevationSampler(dem, data, self.crs) # convert x and y coord of each DEM pixel to the CRS # we are not interpolating the DEM, just taking the points corresponding to each pixel @@ -554,25 +626,15 @@ def toggle(flag): color_on=self.activespace_color ) - def _annotation_z_values(self, linestring: LineString, offset_m: float = 2.0) -> np.ndarray: - """Per-vertex z for airborne tracks that may pass below the local DEM surface.""" - z_vals = [] - for coord in np.array(linestring.coords): - x, y = coord[0], coord[1] - z = vertex_z_from_coord(coord) - lon, lat = self._to_wgs84.transform(x, y) - dem_z = safe_dem_elevation(self.dem, lon, lat) - if z <= 0 or z < dem_z: - z_vals.append(max(dem_z, 0.0) + offset_m) - else: - z_vals.append(z) - return np.asarray(z_vals) - def _annotation_polyline(self, linestring: LineString) -> pv.PolyData: """Build a 3D polyline for one annotation segment.""" - if is_surface_track(np.array(linestring.coords)): + coords = np.array(linestring.coords) + if is_surface_track(coords): return self._flat_sea_surface_polyline(linestring, self.sea_surface_offset_m) - return create_polyline_3d(linestring, z=self._annotation_z_values(linestring)) + return create_polyline_3d( + linestring, + z=annotation_z_profile(linestring, self._dem_sampler), + ) def _sea_surface_polyline(self, linestring: LineString) -> pv.PolyData: """Build a 3D polyline that follows the local water/ground surface.""" diff --git a/tests/scripts/test_viz.py b/tests/scripts/test_viz.py index aa20580..02cd8db 100644 --- a/tests/scripts/test_viz.py +++ b/tests/scripts/test_viz.py @@ -7,6 +7,7 @@ from nps_active_space.scripts.viz import ( Visualizer, + annotation_z_profile, create_polyline_3d, densify_linestring, format_annotation_summary, @@ -19,6 +20,7 @@ resolve_viz_plot_flags, sea_surface_z_profile, track_points_to_linestring, + utm_orientation_axes_kwargs, vertex_z_from_coord, ) @@ -171,14 +173,36 @@ def test_includes_track_and_surface_counts(self): assert "1 elevated" in summary +class TestAnnotationZProfile: + def test_airborne_uses_stored_z_without_dem(self): + line = LineString([(0, 0, 1500.0), (1, 1, 2000.0), (2, 2, 1800.0)]) + + class _FailSampler: + def sample_utm_many(self, x, y): + raise AssertionError("DEM sampling should be skipped for airborne tracks") + + z = annotation_z_profile(line, _FailSampler()) + assert z.tolist() == [1500.0, 2000.0, 1800.0] + + def test_surface_vertices_clamped_to_dem(self): + line = LineString([(0, 0, 0.0), (1, 1, 250.0)]) + + class _FakeSampler: + def sample_utm_many(self, x, y): + return np.full(len(x), 100.0) + + z = annotation_z_profile(line, _FakeSampler(), offset_m=2.0) + assert z[0] == pytest.approx(102.0) + assert z[1] == pytest.approx(250.0) + + class TestOrientationWidgets: - def test_compass_viewport_sits_right_of_axes(self): - ax = Visualizer.axes_viewport - compass = Visualizer.compass_viewport - assert ax[0] == 0.0 - assert compass[0] == pytest.approx(ax[2]) - assert compass[2] > compass[0] - assert compass[3] == pytest.approx(ax[3]) + def test_utm_axes_use_east_north_up_labels(self): + kwargs = utm_orientation_axes_kwargs() + assert kwargs["xlabel"] == "E" + assert kwargs["ylabel"] == "N" + assert kwargs["zlabel"] == "Z" + assert kwargs["viewport"] == (0, 0, 0.2, 0.2) class TestFlatSeaSurfacePolyline: From d706b07cf6d6998ca23b44ee09197aebb7434f24 Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 12:38:53 -0800 Subject: [PATCH 11/17] line width thinner again --- nps_active_space/scripts/viz.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index 81e5379..f971056 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -431,7 +431,7 @@ def _status(self, message: str) -> None: self.logger.info(message) def _add_track_line( - self, polyline: pv.PolyData, *, color: str, line_width: int = 4 + self, polyline: pv.PolyData, *, color: str, line_width: int = 2 ): """Add a track/annotation polyline that reads clearly over the DEM.""" return self.plotter.add_mesh( @@ -439,11 +439,11 @@ def _add_track_line( color=color, line_width=line_width, render_lines_as_tubes=True, - point_size=4, + point_size=2, ) def _add_annotation_lines( - self, polylines: list[pv.PolyData], *, color: str, line_width: int = 3 + self, polylines: list[pv.PolyData], *, color: str, line_width: int = 2 ): """Add many annotation segments as one mesh (much faster than per-segment tubes).""" polylines = [p for p in polylines if p.n_points >= 2] @@ -454,6 +454,7 @@ def _add_annotation_lines( mesh, color=color, line_width=line_width, + point_size=2, render_lines_as_tubes=False, ) @@ -855,7 +856,7 @@ def plot_vessel_tracks( if line.is_empty or line.length == 0: continue polyline = self._sea_surface_polyline(line) - actor = self._add_track_line(polyline, color=self.vessel_track_color, line_width=5) + actor = self._add_track_line(polyline, color=self.vessel_track_color) actors.append(actor) if not actors: From b678406bfae65489a8b9a0cb76413c8e634c490d Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 12:46:28 -0800 Subject: [PATCH 12/17] resturcutre into multiple files --- nps_active_space/scripts/viz.py | 1019 ++------------------------- nps_active_space/viz/__init__.py | 50 ++ nps_active_space/viz/annotations.py | 32 + nps_active_space/viz/cli.py | 194 +++++ nps_active_space/viz/elevation.py | 121 ++++ nps_active_space/viz/geometry.py | 104 +++ nps_active_space/viz/markers.py | 17 + nps_active_space/viz/visualizer.py | 576 +++++++++++++++ tests/{scripts => viz}/test_viz.py | 14 +- 9 files changed, 1148 insertions(+), 979 deletions(-) create mode 100644 nps_active_space/viz/__init__.py create mode 100644 nps_active_space/viz/annotations.py create mode 100644 nps_active_space/viz/cli.py create mode 100644 nps_active_space/viz/elevation.py create mode 100644 nps_active_space/viz/geometry.py create mode 100644 nps_active_space/viz/markers.py create mode 100644 nps_active_space/viz/visualizer.py rename tests/{scripts => viz}/test_viz.py (95%) diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index f971056..fca8044 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -1,979 +1,46 @@ -from __future__ import annotations - -import argparse -import glob -import os -import sys -from datetime import datetime -from pathlib import Path - -import geopandas as gpd -import numpy as np -import pandas as pd -import pyproj -import pyvista as pv -import rasterio -from tqdm import tqdm -from shapely.geometry import LineString, MultiLineString, MultiPolygon, Point, Polygon, box -from nps_active_space.utils.ais import query_ais_mxak -from nps_active_space.utils.helpers import ( - get_deployment, - get_elevation, - get_logger, - load_annotations, - load_DEM, - load_layered_activespace, - load_studyarea, +"""CLI entry: python -m nps_active_space.scripts.viz or scripts/viz.py path.""" + +from nps_active_space.viz import ( + DemElevationSampler, + Visualizer, + annotation_z_profile, + create_polyline_3d, + densify_linestring, + format_annotation_summary, + is_surface_track, + iter_plot_linestrings, + parse_deployment, + parse_existing_file, + parse_iso_date, + parse_max_tracks, + resolve_viz_plot_flags, + sea_surface_z_profile, + track_points_to_linestring, + utm_orientation_axes_kwargs, + vertex_z_from_coord, ) -from nps_active_space.scripts.run_audible_transits import AudibleTransits -from nps_active_space.utils.models import Annotations, Tracks -from nps_active_space.utils.computation import NMSIM_bbox_utm -import nps_active_space.utils.config as cfg - -# CLI argument validation ===================================== - -def parse_deployment(name: str) -> tuple[str, str, int]: - """Parse deployment name like DENATRLA2024 into unit, site, and year.""" - if len(name) < 9: - raise argparse.ArgumentTypeError( - f"deployment must be at least 9 characters " - f"(4-char unit + site code + 4-digit year), got {len(name)}: {name!r}" - ) - unit, site, year_str = name[:4], name[4:-4], name[-4:] - if not site: - raise argparse.ArgumentTypeError( - f"deployment must include a site code between unit and year, got {name!r}" - ) - if len(year_str) != 4 or not year_str.isdigit(): - raise argparse.ArgumentTypeError( - f"deployment year must be 4 digits, got {year_str!r} in {name!r}" - ) - return unit, site, int(year_str) - - -def parse_iso_date(value: str, *, arg_name: str) -> str: - """Validate YYYY-MM-DD date string.""" - try: - datetime.strptime(value, "%Y-%m-%d") - except ValueError: - raise argparse.ArgumentTypeError( - f"{arg_name}: expected YYYY-MM-DD date, got {value!r}" - ) from None - return value - - -def parse_existing_file(path: str, *, arg_name: str) -> str: - """Validate that a CLI path refers to an existing file.""" - file_path = Path(path) - if not file_path.is_file(): - raise argparse.ArgumentTypeError(f"{arg_name}: file not found: {path}") - return path - - -def parse_max_tracks(value: str) -> int: - """Validate positive integer for --max-tracks.""" - try: - n = int(value) - except ValueError: - raise argparse.ArgumentTypeError( - f"--max-tracks must be a positive integer, got {value!r}" - ) from None - if n <= 0: - raise argparse.ArgumentTypeError( - f"--max-tracks must be a positive integer, got {n}" - ) - return n - - -def resolve_viz_plot_flags( - *, - active_space: bool = False, - annotations: bool = False, - audible_transits: bool = False, - vessels: bool = False, - plot_all: bool = False, - annotation_file: str | None = None, - transits_pkl: str | None = None, -) -> tuple[bool, bool, bool, bool]: - """Map CLI flags to Visualizer layer toggles.""" - do_active = active_space or plot_all - do_annotations = annotations or plot_all or annotation_file is not None - do_transits = audible_transits or plot_all or transits_pkl is not None - do_vessels = vessels - return do_active, do_annotations, do_transits, do_vessels - - -def utm_orientation_axes_kwargs() -> dict: - """Orientation-marker kwargs for UTM scenes (+X east, +Y north). - - PyVista/VTK allow only one vtkOrientationMarkerWidget per renderer, so axes - and north-arrow widgets cannot coexist as separate corner widgets. - """ - return { - "interactive": False, - "line_width": 2, - "xlabel": "E", - "ylabel": "N", - "zlabel": "Z", - "viewport": (0, 0, 0.2, 0.2), - } - - -# helper functions ============================================ - -def polygon_to_mesh(polygon: Polygon, z: float) -> pv.PolyData: - exterior = np.array(polygon.exterior.coords) - points = np.c_[exterior[:, 0], exterior[:, 1], np.full(len(exterior), z)] - poly = pv.PolyData(points) - poly.faces = np.hstack([[len(points)], np.arange(len(points))]) - return poly.triangulate() - - -def active_to_polys(active: gpd.GeoDataFrame) -> list[Polygon]: - geometry = active.geometry.iloc[0] - if isinstance(geometry, Polygon): - polys = [geometry] - elif isinstance(geometry, MultiPolygon): - polys = geometry.geoms - return polys - - -def active_to_linestrings(active: gpd.GeoDataFrame) -> list[LineString]: - geometry = active.boundary.iloc[0] - if isinstance(geometry, MultiLineString): - linestrings = geometry.geoms - elif isinstance(geometry, LineString): - linestrings = [geometry] - return linestrings - - -def vertex_z_from_coord(coord: tuple | np.ndarray) -> float: - """Return vertex elevation; missing or NaN z is treated as sea level (0 m).""" - if len(coord) < 3: - return 0.0 - z = float(coord[2]) - return 0.0 if not np.isfinite(z) else z - - -def is_surface_track(coords: np.ndarray) -> bool: - """True when every vertex is at or below sea level (vessel / surface transit).""" - return all(vertex_z_from_coord(c) <= 0.0 for c in coords) - - -def is_airborne_track(coords: np.ndarray) -> bool: - """True when every vertex has positive stored elevation (aircraft spline z).""" - return all(vertex_z_from_coord(c) > 0.0 for c in coords) - - -def stored_vertex_z(linestring: LineString) -> np.ndarray: - """Elevation from annotation geometry coordinates.""" - return np.array([vertex_z_from_coord(c) for c in linestring.coords]) - - -class DemElevationSampler: - """Sample elevations from an in-memory DEM band (no per-point raster I/O).""" - - def __init__(self, dem, band: np.ndarray, plot_crs: str) -> None: - self.dem = dem - self.band = np.asarray(band, dtype=float).copy() - if dem.nodata is not None: - self.band[self.band == dem.nodata] = 0.0 - self.band[self.band > 9000] = 0.0 - self._to_dem = pyproj.Transformer.from_crs(plot_crs, dem.crs, always_xy=True) - - def sample_utm_many(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: - x = np.asarray(x, dtype=float) - y = np.asarray(y, dtype=float) - dem_x, dem_y = self._to_dem.transform(x, y) - rows, cols = self.dem.index(dem_x, dem_y) - rows = np.atleast_1d(rows).astype(int) - cols = np.atleast_1d(cols).astype(int) - in_bounds = ( - (rows >= 0) - & (rows < self.band.shape[0]) - & (cols >= 0) - & (cols < self.band.shape[1]) - ) - elev = np.zeros(rows.shape[0], dtype=float) - if in_bounds.any(): - samples = self.band[rows[in_bounds], cols[in_bounds]] - elev[in_bounds] = np.where(np.isfinite(samples), samples, 0.0) - return elev - - -def annotation_z_profile( - linestring: LineString, - dem_sampler: DemElevationSampler | None, - *, - offset_m: float = 2.0, -) -> np.ndarray: - """Per-vertex z for annotation segments. - - Aircraft splines already carry MSL elevation in geometry; only vertices at or - below sea level need a local DEM clamp (vessels use the flat sea-surface path). - """ - z_stored = stored_vertex_z(linestring) - if np.all(z_stored > 0): - return z_stored - if dem_sampler is None: - return z_stored - coords = np.array(linestring.coords) - need_dem = z_stored <= 0 - if not need_dem.any(): - return z_stored - dem_z = dem_sampler.sample_utm_many(coords[need_dem, 0], coords[need_dem, 1]) - z_out = z_stored.copy() - for j, i in enumerate(np.where(need_dem)[0]): - z = z_stored[i] - dz = float(dem_z[j]) - if z <= 0 or z < dz: - z_out[i] = max(dz, 0.0) + offset_m - return z_out - - -def safe_dem_elevation(dem, lon: float, lat: float) -> float: - """Sample DEM elevation, returning 0 m when the point is off-raster or nodata.""" - try: - elev = float(get_elevation(dem, lon, lat)) - except (IndexError, ValueError, TypeError): - return 0.0 - return elev if np.isfinite(elev) else 0.0 - - -def _point_as_line(point: Point, span_m: float = 10.0) -> LineString: - """Convert a single-point geometry into a short renderable segment.""" - x, y = point.x, point.y - return LineString([(x, y), (x + span_m, y)]) - - -def iter_plot_linestrings(geometry) -> list[LineString]: - """Normalize annotation geometries to LineStrings PyVista can draw.""" - if geometry is None or geometry.is_empty: - return [] - if isinstance(geometry, Point): - return [_point_as_line(geometry)] - if isinstance(geometry, LineString): - if len(geometry.coords) < 2: - return [_point_as_line(Point(geometry.coords[0]))] - return [geometry] - if isinstance(geometry, MultiLineString): - lines: list[LineString] = [] - for part in geometry.geoms: - lines.extend(iter_plot_linestrings(part)) - return lines - return [] - - -def format_annotation_summary(annotations: gpd.GeoDataFrame) -> str: - """One-line summary of a loaded annotations table for logging.""" - if annotations.empty: - return "0 segments" - n_surface = 0 - n_elevated = 0 - for geom in annotations.geometry: - if geom is None or geom.is_empty: - continue - for line in iter_plot_linestrings(geom): - if is_surface_track(np.array(line.coords)): - n_surface += 1 - else: - n_elevated += 1 - geom_types = ", ".join( - f"{k}={v}" for k, v in annotations.geometry.geom_type.value_counts().items() - ) - n_audible = int(annotations["audible"].sum()) if "audible" in annotations.columns else 0 - return ( - f"{len(annotations)} segments, {annotations['_id'].nunique()} tracks, " - f"{n_audible} audible rows, {n_surface} sea-surface / {n_elevated} elevated line(s), " - f"CRS={annotations.crs}, types: {geom_types}" - ) - - -def densify_linestring(linestring: LineString, step_m: float) -> LineString: - """Insert vertices along a line so elevation can be sampled between sparse AIS points.""" - if linestring.is_empty or linestring.length <= step_m: - return linestring - distances = np.arange(0, linestring.length, step_m) - if distances[-1] < linestring.length: - distances = np.append(distances, linestring.length) - return LineString([linestring.interpolate(float(d)) for d in distances]) - - -def sea_surface_z_profile( - linestring: LineString, - dem, - crs: str, - *, - offset_m: float = 5.0, - densify_step_m: float = 100.0, -) -> tuple[LineString, np.ndarray]: - """Return a densified line and z values hugging the DEM water surface (≈0 m).""" - line = densify_linestring(linestring, densify_step_m) - to_wgs84 = pyproj.Transformer.from_crs(crs, "epsg:4326", always_xy=True) - z_vals = [] - for x, y in np.array(line.coords)[:, :2]: - lon, lat = to_wgs84.transform(x, y) - dem_z = safe_dem_elevation(dem, lon, lat) - # NMSIM DEM uses 0 m over water; keep tracks slightly above the mesh to avoid z-fighting. - z_vals.append(max(dem_z, 0.0) + offset_m) - return line, np.asarray(z_vals) - - -def create_polyline_3d( - linestring: LineString, - z: float | np.ndarray | None = None, -) -> pv.PolyData: - coords = np.array(linestring.coords) - xy = coords[:, :2] - if z is not None: - z_arr = np.atleast_1d(np.asarray(z, dtype=float)) - if z_arr.size != coords.shape[0]: - raise ValueError( - f"z must be scalar or have one value per vertex ({coords.shape[0]}), got {z_arr.size}" - ) - coords = np.column_stack((xy, z_arr)) - elif coords.shape[1] == 2: - coords = np.column_stack((xy, np.zeros(coords.shape[0]))) - assert coords.shape[1] == 3 - coords[:, 2] = np.clip(np.nan_to_num(coords[:, 2], nan=0.0), 0, 10000) - n_points = coords.shape[0] - lines = np.hstack([[n_points], np.arange(n_points)]) - poly = pv.PolyData(coords) - poly.lines = lines - return poly - - -def track_points_to_linestring(track_points: gpd.GeoDataFrame) -> LineString: - """Connect ordered track points into a 2D line in the plot CRS.""" - coords = [(geom.x, geom.y) for geom in track_points.geometry] - if len(coords) < 2: - return LineString(coords) - return LineString(coords) - - - - -# main class ===================================================== - -class Visualizer(): - # color config - activespace_color = "orange" - mic_color = "white" - audible_annotation_color = "deepskyblue" - inaudible_annotation_color = "red" - audible_transits_color = "purple" - vessel_track_color = "cyan" - z_scale_toggle_color = "black" # button toggling z scale - sea_surface_offset_m = 5.0 - sea_surface_densify_step_m = 100.0 - - def __init__( - self, - unit: str, - site: str, - year: int, - env: str, - do_active: bool = False, - gain: float | None = None, - do_annots: bool = False, - do_transits: bool = False, - do_vessels: bool = False, - annotation_file: str | None = None, - audible_transits_pkl: str | None = None, - vessel_start_date: str | None = None, - vessel_end_date: str | None = None, - terraced: bool = False, - fill_layers: bool = False, - max_tracks: int = 1000, - ) -> None: - # class metadata - self.unit = unit - self.site = site - self.year = year - self.deployment = f"{unit}{site}{year}" - cfg.initialize(env) - self.project_dir = cfg.read("project", "dir") - self.fill_layers = fill_layers - self.max_tracks = max_tracks - self.logger = get_logger("VIZ", verbose=True) - - # study area and crs - self.study_area = load_studyarea(self.project_dir, self.unit, self.site, self.year) - self.crs = NMSIM_bbox_utm(self.study_area) - self.study_area = self.study_area.to_crs(self.crs) - self._to_wgs84 = pyproj.Transformer.from_crs(self.crs, "epsg:4326", always_xy=True) - - # plot each element - self.plotter = pv.Plotter() - self.plot_dem() - self.plot_mic() - if do_active: - self.plot_activespace(terraced, gain) - if do_annots: - self.plot_annotations(annotation_file) - if do_transits: - self.plot_audible_transits(audible_transits_pkl) - if do_vessels: - self.plot_vessel_tracks(vessel_start_date, vessel_end_date) - - # configure plot parameters and display - self.plotter.enable_terrain_style() - self.setup_z_scale() - self.plotter.add_title(f"{unit}{site}{year}", font_size=12) - self.setup_orientation_widgets() - self.plotter.reset_camera() - self.plotter.camera.elevation = 30 - self.plotter.show() - - def _status(self, message: str) -> None: - """Print and log a user-facing status line (visible when run via python).""" - print(message, flush=True) - self.logger.info(message) - - def _add_track_line( - self, polyline: pv.PolyData, *, color: str, line_width: int = 2 - ): - """Add a track/annotation polyline that reads clearly over the DEM.""" - return self.plotter.add_mesh( - polyline, - color=color, - line_width=line_width, - render_lines_as_tubes=True, - point_size=2, - ) - - def _add_annotation_lines( - self, polylines: list[pv.PolyData], *, color: str, line_width: int = 2 - ): - """Add many annotation segments as one mesh (much faster than per-segment tubes).""" - polylines = [p for p in polylines if p.n_points >= 2] - if not polylines: - return None - mesh = polylines[0] if len(polylines) == 1 else pv.merge(polylines) - return self.plotter.add_mesh( - mesh, - color=color, - line_width=line_width, - point_size=2, - render_lines_as_tubes=False, - ) - - @staticmethod - def _flat_sea_surface_polyline( - linestring: LineString, offset_m: float - ) -> pv.PolyData: - """Sea-level track at a constant offset — no DEM resampling.""" - z = np.full(len(linestring.coords), offset_m) - return create_polyline_3d(linestring, z=z) - - def plot_dem(self, show_scalar_bar: bool = False) -> None: - # load DEM - dem = load_DEM(self.project_dir, self.unit, self.site) - self.dem = dem - data = dem.read(1) - if dem.nodata is not None: - data[data == dem.nodata] = 0 - data[data > 9000] = 0 # higher than any elevation on earth, should be nodata - self._dem_sampler = DemElevationSampler(dem, data, self.crs) - - # convert x and y coord of each DEM pixel to the CRS - # we are not interpolating the DEM, just taking the points corresponding to each pixel - # and putting them in the right place in the target CRS - x = np.arange(dem.shape[1]) - y = np.arange(dem.shape[0]) - x, y = np.meshgrid(x, y) - x_coords, y_coords = rasterio.transform.xy(dem.transform, y, x, offset="center") - x_coords = x_coords.reshape(data.shape) - y_coords = y_coords.reshape(data.shape) - transformer = pyproj.Transformer.from_crs(dem.crs, self.crs, always_xy=True) - x_coords, y_coords = transformer.transform(x_coords, y_coords) - - # create structured grid mesh from DEM - mesh = pv.StructuredGrid() - mesh.points = np.c_[x_coords.flatten(), y_coords.flatten(), data.flatten()] - mesh.dimensions = (dem.shape[1], dem.shape[0], 1) - mesh["elevation"] = data.flatten() - - # add mesh to plot - self.plotter.add_mesh(mesh, scalars="elevation", cmap="gist_earth", show_scalar_bar=show_scalar_bar) - - def plot_point(self, x: float, y: float, z: float, color: str = "white") -> None: - point = pv.PolyData(np.array([[x, y, z]])) - self.plotter.add_mesh(point, color=color, point_size=10, render_points_as_spheres=True) - - def plot_mic(self) -> None: - mic = get_deployment(self.project_dir, self.unit, self.site, self.year) - mic = mic.to_crs(self.crs) - self.plot_point(mic.x, mic.y, mic.z, self.mic_color) - - def plot_activespace(self, terraced: bool = False, gain: float | None = None) -> None: - csv_3d_fits = f"{self.project_dir}/fits.csv" - fit_results = pd.read_csv(csv_3d_fits, index_col="Designator") - - if gain is not None: - print(f"Using gain {gain}dB") - elif self.deployment in fit_results.index: - gain = float(fit_results.loc[self.deployment, "1/3rd Octave Gain (F1)"]) - else: - print(f"No fitted active space gain found in {csv_3d_fits}, skipping active space.") - return - - # load activespace and plot the version specified by the user - active_3d = load_layered_activespace(self.project_dir, self.unit, self.site, self.year, - gain, self.crs) - if terraced: - self.plot_terraced_activespace(active_3d) - else: - self.plot_contoured_activespace(active_3d) - - def plot_contoured_activespace(self, active_3d) -> None: - layer_checkboxes = [] - layer_callbacks = [] - for i, (active_z, active) in enumerate(active_3d.activespaces.items()): - if not active.empty: - checkbox, toggle_cb = self.plot_active_layer(active, active_z, i=i) - layer_checkboxes.append(checkbox) - layer_callbacks.append(toggle_cb) - - # make checkbox to toggle all active space layers - def toggle_all_actives(flag): - for box, toggle_cb in zip(layer_checkboxes, layer_callbacks): - box.GetRepresentation().SetState(int(flag)) - toggle_cb(flag) - self.plotter.render() - - self.plotter.add_checkbox_button_widget( - callback=toggle_all_actives, - value=True, - position=(10, 5), - size=35, - color_on=self.activespace_color - ) - - def plot_active_layer(self, active_layer, elevation: float, i: int = 0): - poly_actor = None - if self.fill_layers: - meshes = [] - for poly in active_to_polys(active_layer): - meshes.append(polygon_to_mesh(poly, elevation)) - poly_data = pv.PolyData().merge(meshes) - poly_actor = self.plotter.add_mesh(poly_data, color=self.activespace_color, opacity=0.5) - - line_actors = [] - for line in active_to_linestrings(active_layer): - polyline = create_polyline_3d(line, z=elevation) - actor = self.plotter.add_mesh(polyline, color=self.activespace_color, point_size=2, line_width=2) - line_actors.append(actor) - - def toggle(flag): - if poly_actor is not None: - poly_actor.SetVisibility(flag) - for actor in line_actors: - actor.SetVisibility(flag) - - checkbox = self.plotter.add_checkbox_button_widget( - callback=toggle, - value=True, - position=(60 + 40*i,10), - size=25, - color_on=self.activespace_color - ) - - return checkbox, toggle - - def plot_terraced_activespace( - self, active_3d, layer_thickness: int = 300, opacity: float = 1 - ) -> None: - meshes = [] - layers = list(active_3d.activespaces.items()) - for i in range(len(layers)): - active_z, active = layers[i] - if active.empty: - continue - - # add the vertical "walls" between the layers - for poly in active_to_polys(active): - hole_polys = [Polygon(hole) for hole in poly.interiors] - for p in [poly] + hole_polys: - mesh = polygon_to_mesh(p, active_z - 0.5 * layer_thickness) - # Extrude the polygon upward by e.g. layer thickness - extruded = mesh.extrude([0, 0, layer_thickness], capping=False) - meshes.append(extruded) - - # add the flat part between the layers - the "floor" of this layer - if i == 0: - floor = active - else: - prev_active = layers[i-1][1].to_crs(self.crs) - sym_diff = active.union_all().symmetric_difference(prev_active.union_all()) - floor = gpd.GeoDataFrame(geometry=[sym_diff], crs=self.crs) - - for poly in active_to_polys(floor): - mesh = polygon_to_mesh(poly, active_z - 0.5*layer_thickness) - meshes.append(mesh) - - # combine all meshes and add to plot - stacked = pv.MultiBlock(meshes).combine() - actor = self.plotter.add_mesh(stacked, color=self.activespace_color, opacity=opacity) - - def toggle(flag): - if actor is not None: - actor.SetVisibility(flag) - self.plotter.add_checkbox_button_widget( - callback=toggle, - value=True, - position=(10, 5), - size=35, - color_on=self.activespace_color - ) - - def _annotation_polyline(self, linestring: LineString) -> pv.PolyData: - """Build a 3D polyline for one annotation segment.""" - coords = np.array(linestring.coords) - if is_surface_track(coords): - return self._flat_sea_surface_polyline(linestring, self.sea_surface_offset_m) - return create_polyline_3d( - linestring, - z=annotation_z_profile(linestring, self._dem_sampler), - ) - - def _sea_surface_polyline(self, linestring: LineString) -> pv.PolyData: - """Build a 3D polyline that follows the local water/ground surface.""" - if is_surface_track(np.array(linestring.coords)): - return self._flat_sea_surface_polyline(linestring, self.sea_surface_offset_m) - line, z_vals = sea_surface_z_profile( - linestring, - self.dem, - self.crs, - offset_m=self.sea_surface_offset_m, - densify_step_m=self.sea_surface_densify_step_m, - ) - return create_polyline_3d(line, z=z_vals) - - def plot_annotations(self, annotation_file: str | None = None) -> None: - if annotation_file is None: - self._status( - f"Loading annotations from project dir for {self.deployment} (valid only)" - ) - annotations = load_annotations( - self.project_dir, self.unit, self.site, self.year, only_valid=True - ) - else: - self._status(f"Loading annotations from {annotation_file} (valid only)") - annotations = Annotations(annotation_file, only_valid=True) - - self._status(f"Parsed annotations: {format_annotation_summary(annotations)}") - if annotations.empty: - self._status("No annotations found, skipping.") - return - - track_ids = annotations["_id"].drop_duplicates() - if len(track_ids) > self.max_tracks: - self._status( - f"Sampling {self.max_tracks} of {len(track_ids)} tracks " - f"({len(annotations)} segments before sample)" - ) - selected_track_ids = track_ids.sample(self.max_tracks, replace=False, random_state=2) - annotations = annotations[annotations["_id"].isin(selected_track_ids)] - self._status(f"After sample: {format_annotation_summary(annotations)}") - - n_loaded = len(annotations) - self._status(f"Reprojecting annotations to {self.crs} and clipping to study area") - annotations = annotations.to_crs(self.crs) - ann_bounds = annotations.total_bounds - study_bounds = self.study_area.total_bounds - self._status(f"Annotation bounds: {ann_bounds}") - self._status(f"Study area bounds: {study_bounds}") - - n_empty = int(annotations.geometry.is_empty.sum()) - if n_empty: - self._status(f"Dropping {n_empty} empty geometries before clip") - annotations = annotations[~annotations.geometry.is_empty] - annotations = annotations.clip(box(*study_bounds)) - annotations = annotations[~annotations.geometry.is_empty].explode(ignore_index=True) - self._status( - f"After clip: {len(annotations)} segments " - f"({n_loaded - len(annotations)} removed outside study area)" - ) - - audible_actors = [] - inaudible_actors = [] - audible_polylines: list[pv.PolyData] = [] - inaudible_polylines: list[pv.PolyData] = [] - n_plotted = 0 - n_skipped = 0 - for _, annot in tqdm( - annotations.iterrows(), - total=len(annotations), - desc="Building annotation lines", - unit="segment", - ): - for line in iter_plot_linestrings(annot["geometry"]): - polyline = self._annotation_polyline(line) - if polyline.n_points < 2: - n_skipped += 1 - continue - n_plotted += 1 - if annot["audible"]: - audible_polylines.append(polyline) - else: - inaudible_polylines.append(polyline) - - audible_actor = self._add_annotation_lines( - audible_polylines, color=self.audible_annotation_color - ) - if audible_actor is not None: - audible_actors.append(audible_actor) - inaudible_actor = self._add_annotation_lines( - inaudible_polylines, color=self.inaudible_annotation_color - ) - if inaudible_actor is not None: - inaudible_actors.append(inaudible_actor) - - if n_plotted == 0: - self._status( - f"No annotation segments plotted ({n_loaded} loaded; " - f"{n_skipped} degenerate lines skipped). " - "Check CRS and study-area overlap." - ) - return - self._status( - f"Plotted {n_plotted} annotation line(s) in " - f"{len(audible_actors) + len(inaudible_actors)} mesh(es) " - f"({len(audible_polylines)} audible, {len(inaudible_polylines)} inaudible; " - f"{n_skipped} skipped)" - ) - - def toggle_audible(flag): - for actor in audible_actors: - actor.SetVisibility(flag) - - def toggle_inaudible(flag): - for actor in inaudible_actors: - actor.SetVisibility(flag) - - self.plotter.add_checkbox_button_widget( - callback=toggle_audible, - value=True, - position=(10,100), - size=25, - color_on="deepskyblue" - ) - self.plotter.add_checkbox_button_widget( - callback=toggle_inaudible, - value=True, - position=(10,60), - size=25, - color_on="red" - ) - - def plot_audible_transits(self, audible_transits_pkl: str | None = None) -> None: - if audible_transits_pkl: - print(f"Loading audible transits from {audible_transits_pkl}") - listener = AudibleTransits.from_pickle(audible_transits_pkl) - else: - print("Loading audible transits") - matches = glob.glob(os.path.join( - self.project_dir, self.unit+self.site, "Output_Data", "AUDIBLE_TRANSITS", - f"3D*{self.year}-01-01*Active Space {self.year}*", "AudibleTransits_object.pkl")) - if len(matches) == 0: - print("No audible transits pkl file found found") - return - listener = AudibleTransits.from_pickle(matches[0]) - - tracks = listener.tracks - if tracks.empty: - print("Audible transits is empty, skipping") - return - - # downsample if too many - print(f"{len(tracks)} audible transits") - if len(tracks) > self.max_tracks: - print("Too many, sampling") - tracks = tracks.sample(self.max_tracks, random_state=4) - print(f"Showing {len(tracks)} transits") - - tracks = tracks.to_crs(self.crs) - - # plot - actors = [] - for _, track in tracks.iterrows(): - polyline = create_polyline_3d(track["interp_geometry"]) - actor = self.plotter.add_mesh(polyline, color=self.audible_transits_color, point_size=2, line_width=2) - actors.append(actor) - - def toggle(flag): - for actor in actors: - actor.SetVisibility(flag) - - self.plotter.add_checkbox_button_widget( - callback=toggle, - value=True, - position=(10,180), - size=25, - color_on="purple" - ) - - def plot_vessel_tracks( - self, start_date: str | None = None, end_date: str | None = None - ) -> None: - """Plot MXAK AIS vessel transits at the local sea surface.""" - start_date = start_date or f"{self.year}-01-01" - end_date = end_date or f"{self.year}-12-31" - - try: - ais_path = cfg.read("data", "ais") - except KeyError: - print("No [data] ais path in config, skipping vessel tracks.") - return - - print(f"Querying AIS tracks from {start_date} to {end_date}") - try: - raw_tracks = query_ais_mxak( - ais_path=ais_path, - start_date=start_date, - end_date=end_date, - mask=self.study_area, - ) - except AssertionError as exc: - print(f"No AIS tracks loaded: {exc}") - return - tracks = Tracks(raw_tracks, id_col="event_id", datetime_col="TIME", z_col="altitude") - tracks = tracks.to_crs(self.crs) - - track_ids = tracks["track_id"].drop_duplicates() - print(f"{len(track_ids)} vessel transits ({len(tracks)} points)") - if len(track_ids) > self.max_tracks: - print(f"More than {self.max_tracks}, sampling") - selected = track_ids.sample(self.max_tracks, replace=False, random_state=3) - tracks = tracks[tracks["track_id"].isin(selected)] - print(f"Showing {selected.nunique()} transits") - - actors = [] - for track_id, group in tracks.groupby("track_id", sort=False): - group = group.sort_values("point_dt") - line = track_points_to_linestring(group) - if line.is_empty or line.length == 0: - continue - polyline = self._sea_surface_polyline(line) - actor = self._add_track_line(polyline, color=self.vessel_track_color) - actors.append(actor) - - if not actors: - print("No vessel tracks to plot.") - return - - def toggle(flag): - for actor in actors: - actor.SetVisibility(flag) - - self.plotter.add_checkbox_button_widget( - callback=toggle, - value=True, - position=(10, 220), - size=25, - color_on=self.vessel_track_color, - ) - - def setup_orientation_widgets(self) -> None: - """Bottom-left E/N/Z axes (+Y = north in UTM).""" - self.plotter.add_axes(**utm_orientation_axes_kwargs()) - - def setup_z_scale(self) -> None: - self.plotter.set_scale(1, 1, 2) # easier to make out elevation differences - - # allow toggling though - def toggle_z_scale(flag): - self.plotter.set_scale(1, 1, 2 if flag else 1) - - self.plotter.add_checkbox_button_widget( - callback=toggle_z_scale, - value=True, - position=(10,140), - size=25, - color_on=self.z_scale_toggle_color - ) - - +from nps_active_space.viz.cli import main + +__all__ = [ + "DemElevationSampler", + "Visualizer", + "annotation_z_profile", + "create_polyline_3d", + "densify_linestring", + "format_annotation_summary", + "is_surface_track", + "iter_plot_linestrings", + "main", + "parse_deployment", + "parse_existing_file", + "parse_iso_date", + "parse_max_tracks", + "resolve_viz_plot_flags", + "sea_surface_z_profile", + "track_points_to_linestring", + "utm_orientation_axes_kwargs", + "vertex_z_from_coord", +] if __name__ == "__main__": - parser = argparse.ArgumentParser() - - # required args - # we pass deployment as a single positional arg instead of -u DENA -s TRLA -y 2024, - # because it saves a lot of keystrokes, and this script is often used frequently - parser.add_argument( - "deployment", - type=parse_deployment, - help="Deployment name, e.g. DENATRLA2024", - ) - parser.add_argument("-e", "--environment", default="DENA_streamline", - help="Config environment name, e.g. DENA_streamline. Default is 'DENA_streamline'") - - # common args - parser.add_argument("-s", "--active-space", action="store_true", - help="If included, load and plot the active space.") - parser.add_argument("-g", "--gain", type=float, help="Active space gain, if not the default.") - parser.add_argument("-a", "--annotations", action="store_true", - help="If included, load and plot annotations") - parser.add_argument("-t", "--audible-transits", action="store_true", - help="If included, load and plot audible transits") - parser.add_argument("-v", "--vessels", action="store_true", - help="If included, load and plot MXAK AIS vessel tracks at sea level") - parser.add_argument("--all", action="store_true", - help="Load everything, shorthand for --active-space --annotations --audible-transits") - - # uncommon / special use case args - parser.add_argument( - "-m", - "--max-tracks", - type=parse_max_tracks, - default=500, - help="Maximum number of annotation tracks or audible transits to show.", - ) - parser.add_argument( - "--annotation-file", - type=lambda p: parse_existing_file(p, arg_name="--annotation-file"), - help="Path to .geojson annotations (implies -a).", - ) - parser.add_argument( - "--transits-pkl", - type=lambda p: parse_existing_file(p, arg_name="--transits-pkl"), - help="Path to .pkl audible transits (implies -t).", - ) - parser.add_argument( - "--start-date", - type=lambda d: parse_iso_date(d, arg_name="--start-date"), - help="AIS query start date (YYYY-MM-DD). Default: Jan 1 of deployment year.", - ) - parser.add_argument( - "--end-date", - type=lambda d: parse_iso_date(d, arg_name="--end-date"), - help="AIS query end date (YYYY-MM-DD). Default: Dec 31 of deployment year.", - ) - parser.add_argument("--terraced", action="store_true", - help="If included, render the active space as the terraced surface instead of contours.") - parser.add_argument("--fill-layers", action="store_true", - help="If included, fill the interior of each active space contour polygon.") - - args = parser.parse_args() - - unit, site, year = args.deployment - print(unit, site, year) - - do_active, do_annotations, do_transits, do_vessels = resolve_viz_plot_flags( - active_space=args.active_space, - annotations=args.annotations, - audible_transits=args.audible_transits, - vessels=args.vessels, - plot_all=args.all, - annotation_file=args.annotation_file, - transits_pkl=args.transits_pkl, - ) - - Visualizer(unit, site, year, args.environment, - do_active, args.gain, - do_annotations, do_transits, do_vessels, - args.annotation_file, args.transits_pkl, - args.start_date, args.end_date, - args.terraced, args.fill_layers, args.max_tracks) \ No newline at end of file + main() diff --git a/nps_active_space/viz/__init__.py b/nps_active_space/viz/__init__.py new file mode 100644 index 0000000..dbf658e --- /dev/null +++ b/nps_active_space/viz/__init__.py @@ -0,0 +1,50 @@ +"""3D deployment visualization (PyVista).""" + +from nps_active_space.viz.annotations import format_annotation_summary +from nps_active_space.viz.cli import ( + main, + parse_deployment, + parse_existing_file, + parse_iso_date, + parse_max_tracks, + resolve_viz_plot_flags, +) +from nps_active_space.viz.markers import utm_orientation_axes_kwargs +from nps_active_space.viz.elevation import ( + DemElevationSampler, + annotation_z_profile, + is_airborne_track, + is_surface_track, + sea_surface_z_profile, + vertex_z_from_coord, +) +from nps_active_space.viz.geometry import ( + create_polyline_3d, + densify_linestring, + flat_sea_surface_polyline, + iter_plot_linestrings, + track_points_to_linestring, +) +from nps_active_space.viz.visualizer import Visualizer + +__all__ = [ + "DemElevationSampler", + "Visualizer", + "annotation_z_profile", + "create_polyline_3d", + "densify_linestring", + "format_annotation_summary", + "is_airborne_track", + "is_surface_track", + "iter_plot_linestrings", + "main", + "parse_deployment", + "parse_existing_file", + "parse_iso_date", + "parse_max_tracks", + "resolve_viz_plot_flags", + "sea_surface_z_profile", + "track_points_to_linestring", + "utm_orientation_axes_kwargs", + "vertex_z_from_coord", +] diff --git a/nps_active_space/viz/annotations.py b/nps_active_space/viz/annotations.py new file mode 100644 index 0000000..e6f50de --- /dev/null +++ b/nps_active_space/viz/annotations.py @@ -0,0 +1,32 @@ +from __future__ import annotations + +import geopandas as gpd +import numpy as np + +from nps_active_space.viz.elevation import is_surface_track +from nps_active_space.viz.geometry import iter_plot_linestrings + + +def format_annotation_summary(annotations: gpd.GeoDataFrame) -> str: + """One-line summary of a loaded annotations table for logging.""" + if annotations.empty: + return "0 segments" + n_surface = 0 + n_elevated = 0 + for geom in annotations.geometry: + if geom is None or geom.is_empty: + continue + for line in iter_plot_linestrings(geom): + if is_surface_track(np.array(line.coords)): + n_surface += 1 + else: + n_elevated += 1 + geom_types = ", ".join( + f"{k}={v}" for k, v in annotations.geometry.geom_type.value_counts().items() + ) + n_audible = int(annotations["audible"].sum()) if "audible" in annotations.columns else 0 + return ( + f"{len(annotations)} segments, {annotations['_id'].nunique()} tracks, " + f"{n_audible} audible rows, {n_surface} sea-surface / {n_elevated} elevated line(s), " + f"CRS={annotations.crs}, types: {geom_types}" + ) diff --git a/nps_active_space/viz/cli.py b/nps_active_space/viz/cli.py new file mode 100644 index 0000000..bc5c773 --- /dev/null +++ b/nps_active_space/viz/cli.py @@ -0,0 +1,194 @@ +from __future__ import annotations + +import argparse +from datetime import datetime +from pathlib import Path + +from nps_active_space.viz.visualizer import Visualizer + + +def parse_deployment(name: str) -> tuple[str, str, int]: + """Parse deployment name like DENATRLA2024 into unit, site, and year.""" + if len(name) < 9: + raise argparse.ArgumentTypeError( + f"deployment must be at least 9 characters " + f"(4-char unit + site code + 4-digit year), got {len(name)}: {name!r}" + ) + unit, site, year_str = name[:4], name[4:-4], name[-4:] + if not site: + raise argparse.ArgumentTypeError( + f"deployment must include a site code between unit and year, got {name!r}" + ) + if len(year_str) != 4 or not year_str.isdigit(): + raise argparse.ArgumentTypeError( + f"deployment year must be 4 digits, got {year_str!r} in {name!r}" + ) + return unit, site, int(year_str) + + +def parse_iso_date(value: str, *, arg_name: str) -> str: + """Validate YYYY-MM-DD date string.""" + try: + datetime.strptime(value, "%Y-%m-%d") + except ValueError: + raise argparse.ArgumentTypeError( + f"{arg_name}: expected YYYY-MM-DD date, got {value!r}" + ) from None + return value + + +def parse_existing_file(path: str, *, arg_name: str) -> str: + """Validate that a CLI path refers to an existing file.""" + file_path = Path(path) + if not file_path.is_file(): + raise argparse.ArgumentTypeError(f"{arg_name}: file not found: {path}") + return path + + +def parse_max_tracks(value: str) -> int: + """Validate positive integer for --max-tracks.""" + try: + n = int(value) + except ValueError: + raise argparse.ArgumentTypeError( + f"--max-tracks must be a positive integer, got {value!r}" + ) from None + if n <= 0: + raise argparse.ArgumentTypeError( + f"--max-tracks must be a positive integer, got {n}" + ) + return n + + +def resolve_viz_plot_flags( + *, + active_space: bool = False, + annotations: bool = False, + audible_transits: bool = False, + vessels: bool = False, + plot_all: bool = False, + annotation_file: str | None = None, + transits_pkl: str | None = None, +) -> tuple[bool, bool, bool, bool]: + """Map CLI flags to Visualizer layer toggles.""" + do_active = active_space or plot_all + do_annotations = annotations or plot_all or annotation_file is not None + do_transits = audible_transits or plot_all or transits_pkl is not None + do_vessels = vessels + return do_active, do_annotations, do_transits, do_vessels + + +def main() -> None: + parser = argparse.ArgumentParser() + + parser.add_argument( + "deployment", + type=parse_deployment, + help="Deployment name, e.g. DENATRLA2024", + ) + parser.add_argument( + "-e", + "--environment", + default="DENA_streamline", + help="Config environment name, e.g. DENA_streamline. Default is 'DENA_streamline'", + ) + parser.add_argument( + "-s", + "--active-space", + action="store_true", + help="If included, load and plot the active space.", + ) + parser.add_argument("-g", "--gain", type=float, help="Active space gain, if not the default.") + parser.add_argument( + "-a", + "--annotations", + action="store_true", + help="If included, load and plot annotations", + ) + parser.add_argument( + "-t", + "--audible-transits", + action="store_true", + help="If included, load and plot audible transits", + ) + parser.add_argument( + "-v", + "--vessels", + action="store_true", + help="If included, load and plot MXAK AIS vessel tracks at sea level", + ) + parser.add_argument( + "--all", + action="store_true", + help="Load everything, shorthand for --active-space --annotations --audible-transits", + ) + parser.add_argument( + "-m", + "--max-tracks", + type=parse_max_tracks, + default=500, + help="Maximum number of annotation tracks or audible transits to show.", + ) + parser.add_argument( + "--annotation-file", + type=lambda p: parse_existing_file(p, arg_name="--annotation-file"), + help="Path to .geojson annotations (implies -a).", + ) + parser.add_argument( + "--transits-pkl", + type=lambda p: parse_existing_file(p, arg_name="--transits-pkl"), + help="Path to .pkl audible transits (implies -t).", + ) + parser.add_argument( + "--start-date", + type=lambda d: parse_iso_date(d, arg_name="--start-date"), + help="AIS query start date (YYYY-MM-DD). Default: Jan 1 of deployment year.", + ) + parser.add_argument( + "--end-date", + type=lambda d: parse_iso_date(d, arg_name="--end-date"), + help="AIS query end date (YYYY-MM-DD). Default: Dec 31 of deployment year.", + ) + parser.add_argument( + "--terraced", + action="store_true", + help="If included, render the active space as the terraced surface instead of contours.", + ) + parser.add_argument( + "--fill-layers", + action="store_true", + help="If included, fill the interior of each active space contour polygon.", + ) + + args = parser.parse_args() + unit, site, year = args.deployment + print(unit, site, year) + + do_active, do_annotations, do_transits, do_vessels = resolve_viz_plot_flags( + active_space=args.active_space, + annotations=args.annotations, + audible_transits=args.audible_transits, + vessels=args.vessels, + plot_all=args.all, + annotation_file=args.annotation_file, + transits_pkl=args.transits_pkl, + ) + + Visualizer( + unit, + site, + year, + args.environment, + do_active, + args.gain, + do_annotations, + do_transits, + do_vessels, + args.annotation_file, + args.transits_pkl, + args.start_date, + args.end_date, + args.terraced, + args.fill_layers, + args.max_tracks, + ) diff --git a/nps_active_space/viz/elevation.py b/nps_active_space/viz/elevation.py new file mode 100644 index 0000000..a89d777 --- /dev/null +++ b/nps_active_space/viz/elevation.py @@ -0,0 +1,121 @@ +from __future__ import annotations + +import numpy as np +import pyproj +from shapely.geometry import LineString + +from nps_active_space.utils.helpers import get_elevation +from nps_active_space.viz.geometry import densify_linestring + + +def vertex_z_from_coord(coord: tuple | np.ndarray) -> float: + """Return vertex elevation; missing or NaN z is treated as sea level (0 m).""" + if len(coord) < 3: + return 0.0 + z = float(coord[2]) + return 0.0 if not np.isfinite(z) else z + + +def is_surface_track(coords: np.ndarray) -> bool: + """True when every vertex is at or below sea level (vessel / surface transit).""" + return all(vertex_z_from_coord(c) <= 0.0 for c in coords) + + +def is_airborne_track(coords: np.ndarray) -> bool: + """True when every vertex has positive stored elevation (aircraft spline z).""" + return all(vertex_z_from_coord(c) > 0.0 for c in coords) + + +def stored_vertex_z(linestring: LineString) -> np.ndarray: + """Elevation from annotation geometry coordinates.""" + return np.array([vertex_z_from_coord(c) for c in linestring.coords]) + + +class DemElevationSampler: + """Sample elevations from an in-memory DEM band (no per-point raster I/O).""" + + def __init__(self, dem, band: np.ndarray, plot_crs: str) -> None: + self.dem = dem + self.band = np.asarray(band, dtype=float).copy() + if dem.nodata is not None: + self.band[self.band == dem.nodata] = 0.0 + self.band[self.band > 9000] = 0.0 + self._to_dem = pyproj.Transformer.from_crs(plot_crs, dem.crs, always_xy=True) + + def sample_utm_many(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + dem_x, dem_y = self._to_dem.transform(x, y) + rows, cols = self.dem.index(dem_x, dem_y) + rows = np.atleast_1d(rows).astype(int) + cols = np.atleast_1d(cols).astype(int) + in_bounds = ( + (rows >= 0) + & (rows < self.band.shape[0]) + & (cols >= 0) + & (cols < self.band.shape[1]) + ) + elev = np.zeros(rows.shape[0], dtype=float) + if in_bounds.any(): + samples = self.band[rows[in_bounds], cols[in_bounds]] + elev[in_bounds] = np.where(np.isfinite(samples), samples, 0.0) + return elev + + +def annotation_z_profile( + linestring: LineString, + dem_sampler: DemElevationSampler | None, + *, + offset_m: float = 2.0, +) -> np.ndarray: + """Per-vertex z for annotation segments. + + Aircraft splines already carry MSL elevation in geometry; only vertices at or + below sea level need a local DEM clamp (vessels use the flat sea-surface path). + """ + z_stored = stored_vertex_z(linestring) + if np.all(z_stored > 0): + return z_stored + if dem_sampler is None: + return z_stored + coords = np.array(linestring.coords) + need_dem = z_stored <= 0 + if not need_dem.any(): + return z_stored + dem_z = dem_sampler.sample_utm_many(coords[need_dem, 0], coords[need_dem, 1]) + z_out = z_stored.copy() + for j, i in enumerate(np.where(need_dem)[0]): + z = z_stored[i] + dz = float(dem_z[j]) + if z <= 0 or z < dz: + z_out[i] = max(dz, 0.0) + offset_m + return z_out + + +def safe_dem_elevation(dem, lon: float, lat: float) -> float: + """Sample DEM elevation, returning 0 m when the point is off-raster or nodata.""" + try: + elev = float(get_elevation(dem, lon, lat)) + except (IndexError, ValueError, TypeError): + return 0.0 + return elev if np.isfinite(elev) else 0.0 + + +def sea_surface_z_profile( + linestring: LineString, + dem, + crs: str, + *, + offset_m: float = 5.0, + densify_step_m: float = 100.0, +) -> tuple[LineString, np.ndarray]: + """Return a densified line and z values hugging the DEM water surface (≈0 m).""" + line = densify_linestring(linestring, densify_step_m) + to_wgs84 = pyproj.Transformer.from_crs(crs, "epsg:4326", always_xy=True) + z_vals = [] + for x, y in np.array(line.coords)[:, :2]: + lon, lat = to_wgs84.transform(x, y) + dem_z = safe_dem_elevation(dem, lon, lat) + # NMSIM DEM uses 0 m over water; keep tracks slightly above the mesh to avoid z-fighting. + z_vals.append(max(dem_z, 0.0) + offset_m) + return line, np.asarray(z_vals) diff --git a/nps_active_space/viz/geometry.py b/nps_active_space/viz/geometry.py new file mode 100644 index 0000000..bd066ff --- /dev/null +++ b/nps_active_space/viz/geometry.py @@ -0,0 +1,104 @@ +from __future__ import annotations + +import geopandas as gpd +import numpy as np +import pyvista as pv +from shapely.geometry import LineString, MultiLineString, MultiPolygon, Point, Polygon + + +def polygon_to_mesh(polygon: Polygon, z: float) -> pv.PolyData: + exterior = np.array(polygon.exterior.coords) + points = np.c_[exterior[:, 0], exterior[:, 1], np.full(len(exterior), z)] + poly = pv.PolyData(points) + poly.faces = np.hstack([[len(points)], np.arange(len(points))]) + return poly.triangulate() + + +def active_to_polys(active: gpd.GeoDataFrame) -> list[Polygon]: + geometry = active.geometry.iloc[0] + if isinstance(geometry, Polygon): + polys = [geometry] + elif isinstance(geometry, MultiPolygon): + polys = geometry.geoms + return polys + + +def active_to_linestrings(active: gpd.GeoDataFrame) -> list[LineString]: + geometry = active.boundary.iloc[0] + if isinstance(geometry, MultiLineString): + linestrings = geometry.geoms + elif isinstance(geometry, LineString): + linestrings = [geometry] + return linestrings + + +def _point_as_line(point: Point, span_m: float = 10.0) -> LineString: + """Convert a single-point geometry into a short renderable segment.""" + x, y = point.x, point.y + return LineString([(x, y), (x + span_m, y)]) + + +def iter_plot_linestrings(geometry) -> list[LineString]: + """Normalize annotation geometries to LineStrings PyVista can draw.""" + if geometry is None or geometry.is_empty: + return [] + if isinstance(geometry, Point): + return [_point_as_line(geometry)] + if isinstance(geometry, LineString): + if len(geometry.coords) < 2: + return [_point_as_line(Point(geometry.coords[0]))] + return [geometry] + if isinstance(geometry, MultiLineString): + lines: list[LineString] = [] + for part in geometry.geoms: + lines.extend(iter_plot_linestrings(part)) + return lines + return [] + + +def densify_linestring(linestring: LineString, step_m: float) -> LineString: + """Insert vertices along a line so elevation can be sampled between sparse AIS points.""" + if linestring.is_empty or linestring.length <= step_m: + return linestring + distances = np.arange(0, linestring.length, step_m) + if distances[-1] < linestring.length: + distances = np.append(distances, linestring.length) + return LineString([linestring.interpolate(float(d)) for d in distances]) + + +def create_polyline_3d( + linestring: LineString, + z: float | np.ndarray | None = None, +) -> pv.PolyData: + coords = np.array(linestring.coords) + xy = coords[:, :2] + if z is not None: + z_arr = np.atleast_1d(np.asarray(z, dtype=float)) + if z_arr.size != coords.shape[0]: + raise ValueError( + f"z must be scalar or have one value per vertex ({coords.shape[0]}), got {z_arr.size}" + ) + coords = np.column_stack((xy, z_arr)) + elif coords.shape[1] == 2: + coords = np.column_stack((xy, np.zeros(coords.shape[0]))) + assert coords.shape[1] == 3 + coords[:, 2] = np.clip(np.nan_to_num(coords[:, 2], nan=0.0), 0, 10000) + n_points = coords.shape[0] + lines = np.hstack([[n_points], np.arange(n_points)]) + poly = pv.PolyData(coords) + poly.lines = lines + return poly + + +def flat_sea_surface_polyline(linestring: LineString, offset_m: float) -> pv.PolyData: + """Sea-level track at a constant offset — no DEM resampling.""" + z = np.full(len(linestring.coords), offset_m) + return create_polyline_3d(linestring, z=z) + + +def track_points_to_linestring(track_points: gpd.GeoDataFrame) -> LineString: + """Connect ordered track points into a 2D line in the plot CRS.""" + coords = [(geom.x, geom.y) for geom in track_points.geometry] + if len(coords) < 2: + return LineString(coords) + return LineString(coords) diff --git a/nps_active_space/viz/markers.py b/nps_active_space/viz/markers.py new file mode 100644 index 0000000..93e18e8 --- /dev/null +++ b/nps_active_space/viz/markers.py @@ -0,0 +1,17 @@ +"""Scene orientation markers and viewport chrome.""" + + +def utm_orientation_axes_kwargs() -> dict: + """Orientation-marker kwargs for UTM scenes (+X east, +Y north). + + PyVista/VTK allow only one vtkOrientationMarkerWidget per renderer, so axes + and north-arrow widgets cannot coexist as separate corner widgets. + """ + return { + "interactive": False, + "line_width": 2, + "xlabel": "E", + "ylabel": "N", + "zlabel": "Z", + "viewport": (0, 0, 0.2, 0.2), + } diff --git a/nps_active_space/viz/visualizer.py b/nps_active_space/viz/visualizer.py new file mode 100644 index 0000000..d21ff57 --- /dev/null +++ b/nps_active_space/viz/visualizer.py @@ -0,0 +1,576 @@ +from __future__ import annotations + +import glob +import os + +import geopandas as gpd +import numpy as np +import pandas as pd +import pyproj +import pyvista as pv +import rasterio +from shapely.geometry import LineString, Polygon, box +from tqdm import tqdm + +import nps_active_space.utils.config as cfg +from nps_active_space.scripts.run_audible_transits import AudibleTransits +from nps_active_space.utils.ais import query_ais_mxak +from nps_active_space.utils.computation import NMSIM_bbox_utm +from nps_active_space.utils.helpers import ( + get_deployment, + get_logger, + load_annotations, + load_DEM, + load_layered_activespace, + load_studyarea, +) +from nps_active_space.utils.models import Annotations, Tracks +from nps_active_space.viz.annotations import format_annotation_summary +from nps_active_space.viz.markers import utm_orientation_axes_kwargs +from nps_active_space.viz.elevation import ( + DemElevationSampler, + annotation_z_profile, + is_surface_track, + sea_surface_z_profile, +) +from nps_active_space.viz.geometry import ( + active_to_linestrings, + active_to_polys, + create_polyline_3d, + flat_sea_surface_polyline, + iter_plot_linestrings, + polygon_to_mesh, + track_points_to_linestring, +) + + +class Visualizer: + # color config + activespace_color = "orange" + mic_color = "white" + audible_annotation_color = "deepskyblue" + inaudible_annotation_color = "red" + audible_transits_color = "purple" + vessel_track_color = "cyan" + z_scale_toggle_color = "black" # button toggling z scale + sea_surface_offset_m = 5.0 + sea_surface_densify_step_m = 100.0 + + def __init__( + self, + unit: str, + site: str, + year: int, + env: str, + do_active: bool = False, + gain: float | None = None, + do_annots: bool = False, + do_transits: bool = False, + do_vessels: bool = False, + annotation_file: str | None = None, + audible_transits_pkl: str | None = None, + vessel_start_date: str | None = None, + vessel_end_date: str | None = None, + terraced: bool = False, + fill_layers: bool = False, + max_tracks: int = 1000, + ) -> None: + self.unit = unit + self.site = site + self.year = year + self.deployment = f"{unit}{site}{year}" + cfg.initialize(env) + self.project_dir = cfg.read("project", "dir") + self.fill_layers = fill_layers + self.max_tracks = max_tracks + self.logger = get_logger("VIZ", verbose=True) + + self.study_area = load_studyarea(self.project_dir, self.unit, self.site, self.year) + self.crs = NMSIM_bbox_utm(self.study_area) + self.study_area = self.study_area.to_crs(self.crs) + self._to_wgs84 = pyproj.Transformer.from_crs(self.crs, "epsg:4326", always_xy=True) + + self.plotter = pv.Plotter() + self.plot_dem() + self.plot_mic() + if do_active: + self.plot_activespace(terraced, gain) + if do_annots: + self.plot_annotations(annotation_file) + if do_transits: + self.plot_audible_transits(audible_transits_pkl) + if do_vessels: + self.plot_vessel_tracks(vessel_start_date, vessel_end_date) + + self.plotter.enable_terrain_style() + self.setup_z_scale() + self.plotter.add_title(f"{unit}{site}{year}", font_size=12) + self.setup_orientation_widgets() + self.plotter.reset_camera() + self.plotter.camera.elevation = 30 + self.plotter.show() + + def _status(self, message: str) -> None: + """Print and log a user-facing status line (visible when run via python).""" + print(message, flush=True) + self.logger.info(message) + + def _add_track_line( + self, polyline: pv.PolyData, *, color: str, line_width: int = 2 + ): + """Add a track/annotation polyline that reads clearly over the DEM.""" + return self.plotter.add_mesh( + polyline, + color=color, + line_width=line_width, + render_lines_as_tubes=True, + point_size=2, + ) + + def _add_annotation_lines( + self, polylines: list[pv.PolyData], *, color: str, line_width: int = 2 + ): + """Add many annotation segments as one mesh (much faster than per-segment tubes).""" + polylines = [p for p in polylines if p.n_points >= 2] + if not polylines: + return None + mesh = polylines[0] if len(polylines) == 1 else pv.merge(polylines) + return self.plotter.add_mesh( + mesh, + color=color, + line_width=line_width, + point_size=2, + render_lines_as_tubes=False, + ) + + @staticmethod + def _flat_sea_surface_polyline( + linestring: LineString, offset_m: float + ) -> pv.PolyData: + return flat_sea_surface_polyline(linestring, offset_m) + + def plot_dem(self, show_scalar_bar: bool = False) -> None: + dem = load_DEM(self.project_dir, self.unit, self.site) + self.dem = dem + data = dem.read(1) + if dem.nodata is not None: + data[data == dem.nodata] = 0 + data[data > 9000] = 0 + self._dem_sampler = DemElevationSampler(dem, data, self.crs) + + x = np.arange(dem.shape[1]) + y = np.arange(dem.shape[0]) + x, y = np.meshgrid(x, y) + x_coords, y_coords = rasterio.transform.xy(dem.transform, y, x, offset="center") + x_coords = x_coords.reshape(data.shape) + y_coords = y_coords.reshape(data.shape) + transformer = pyproj.Transformer.from_crs(dem.crs, self.crs, always_xy=True) + x_coords, y_coords = transformer.transform(x_coords, y_coords) + + mesh = pv.StructuredGrid() + mesh.points = np.c_[x_coords.flatten(), y_coords.flatten(), data.flatten()] + mesh.dimensions = (dem.shape[1], dem.shape[0], 1) + mesh["elevation"] = data.flatten() + + self.plotter.add_mesh(mesh, scalars="elevation", cmap="gist_earth", show_scalar_bar=show_scalar_bar) + + def plot_point(self, x: float, y: float, z: float, color: str = "white") -> None: + point = pv.PolyData(np.array([[x, y, z]])) + self.plotter.add_mesh(point, color=color, point_size=10, render_points_as_spheres=True) + + def plot_mic(self) -> None: + mic = get_deployment(self.project_dir, self.unit, self.site, self.year) + mic = mic.to_crs(self.crs) + self.plot_point(mic.x, mic.y, mic.z, self.mic_color) + + def plot_activespace(self, terraced: bool = False, gain: float | None = None) -> None: + csv_3d_fits = f"{self.project_dir}/fits.csv" + fit_results = pd.read_csv(csv_3d_fits, index_col="Designator") + + if gain is not None: + print(f"Using gain {gain}dB") + elif self.deployment in fit_results.index: + gain = float(fit_results.loc[self.deployment, "1/3rd Octave Gain (F1)"]) + else: + print(f"No fitted active space gain found in {csv_3d_fits}, skipping active space.") + return + + active_3d = load_layered_activespace( + self.project_dir, self.unit, self.site, self.year, gain, self.crs + ) + if terraced: + self.plot_terraced_activespace(active_3d) + else: + self.plot_contoured_activespace(active_3d) + + def plot_contoured_activespace(self, active_3d) -> None: + layer_checkboxes = [] + layer_callbacks = [] + for i, (active_z, active) in enumerate(active_3d.activespaces.items()): + if not active.empty: + checkbox, toggle_cb = self.plot_active_layer(active, active_z, i=i) + layer_checkboxes.append(checkbox) + layer_callbacks.append(toggle_cb) + + def toggle_all_actives(flag): + for box, toggle_cb in zip(layer_checkboxes, layer_callbacks): + box.GetRepresentation().SetState(int(flag)) + toggle_cb(flag) + self.plotter.render() + + self.plotter.add_checkbox_button_widget( + callback=toggle_all_actives, + value=True, + position=(10, 5), + size=35, + color_on=self.activespace_color, + ) + + def plot_active_layer(self, active_layer, elevation: float, i: int = 0): + poly_actor = None + if self.fill_layers: + meshes = [] + for poly in active_to_polys(active_layer): + meshes.append(polygon_to_mesh(poly, elevation)) + poly_data = pv.PolyData().merge(meshes) + poly_actor = self.plotter.add_mesh(poly_data, color=self.activespace_color, opacity=0.5) + + line_actors = [] + for line in active_to_linestrings(active_layer): + polyline = create_polyline_3d(line, z=elevation) + actor = self.plotter.add_mesh( + polyline, color=self.activespace_color, point_size=2, line_width=2 + ) + line_actors.append(actor) + + def toggle(flag): + if poly_actor is not None: + poly_actor.SetVisibility(flag) + for actor in line_actors: + actor.SetVisibility(flag) + + checkbox = self.plotter.add_checkbox_button_widget( + callback=toggle, + value=True, + position=(60 + 40 * i, 10), + size=25, + color_on=self.activespace_color, + ) + + return checkbox, toggle + + def plot_terraced_activespace( + self, active_3d, layer_thickness: int = 300, opacity: float = 1 + ) -> None: + meshes = [] + layers = list(active_3d.activespaces.items()) + for i in range(len(layers)): + active_z, active = layers[i] + if active.empty: + continue + + for poly in active_to_polys(active): + hole_polys = [Polygon(hole) for hole in poly.interiors] + for p in [poly] + hole_polys: + mesh = polygon_to_mesh(p, active_z - 0.5 * layer_thickness) + extruded = mesh.extrude([0, 0, layer_thickness], capping=False) + meshes.append(extruded) + + if i == 0: + floor = active + else: + prev_active = layers[i - 1][1].to_crs(self.crs) + sym_diff = active.union_all().symmetric_difference(prev_active.union_all()) + floor = gpd.GeoDataFrame(geometry=[sym_diff], crs=self.crs) + + for poly in active_to_polys(floor): + mesh = polygon_to_mesh(poly, active_z - 0.5 * layer_thickness) + meshes.append(mesh) + + stacked = pv.MultiBlock(meshes).combine() + actor = self.plotter.add_mesh(stacked, color=self.activespace_color, opacity=opacity) + + def toggle(flag): + if actor is not None: + actor.SetVisibility(flag) + + self.plotter.add_checkbox_button_widget( + callback=toggle, + value=True, + position=(10, 5), + size=35, + color_on=self.activespace_color, + ) + + def _annotation_polyline(self, linestring: LineString) -> pv.PolyData: + """Build a 3D polyline for one annotation segment.""" + coords = np.array(linestring.coords) + if is_surface_track(coords): + return self._flat_sea_surface_polyline(linestring, self.sea_surface_offset_m) + return create_polyline_3d( + linestring, + z=annotation_z_profile(linestring, self._dem_sampler), + ) + + def _sea_surface_polyline(self, linestring: LineString) -> pv.PolyData: + """Build a 3D polyline that follows the local water/ground surface.""" + if is_surface_track(np.array(linestring.coords)): + return self._flat_sea_surface_polyline(linestring, self.sea_surface_offset_m) + line, z_vals = sea_surface_z_profile( + linestring, + self.dem, + self.crs, + offset_m=self.sea_surface_offset_m, + densify_step_m=self.sea_surface_densify_step_m, + ) + return create_polyline_3d(line, z=z_vals) + + def plot_annotations(self, annotation_file: str | None = None) -> None: + if annotation_file is None: + self._status( + f"Loading annotations from project dir for {self.deployment} (valid only)" + ) + annotations = load_annotations( + self.project_dir, self.unit, self.site, self.year, only_valid=True + ) + else: + self._status(f"Loading annotations from {annotation_file} (valid only)") + annotations = Annotations(annotation_file, only_valid=True) + + self._status(f"Parsed annotations: {format_annotation_summary(annotations)}") + if annotations.empty: + self._status("No annotations found, skipping.") + return + + track_ids = annotations["_id"].drop_duplicates() + if len(track_ids) > self.max_tracks: + self._status( + f"Sampling {self.max_tracks} of {len(track_ids)} tracks " + f"({len(annotations)} segments before sample)" + ) + selected_track_ids = track_ids.sample(self.max_tracks, replace=False, random_state=2) + annotations = annotations[annotations["_id"].isin(selected_track_ids)] + self._status(f"After sample: {format_annotation_summary(annotations)}") + + n_loaded = len(annotations) + self._status(f"Reprojecting annotations to {self.crs} and clipping to study area") + annotations = annotations.to_crs(self.crs) + ann_bounds = annotations.total_bounds + study_bounds = self.study_area.total_bounds + self._status(f"Annotation bounds: {ann_bounds}") + self._status(f"Study area bounds: {study_bounds}") + + n_empty = int(annotations.geometry.is_empty.sum()) + if n_empty: + self._status(f"Dropping {n_empty} empty geometries before clip") + annotations = annotations[~annotations.geometry.is_empty] + annotations = annotations.clip(box(*study_bounds)) + annotations = annotations[~annotations.geometry.is_empty].explode(ignore_index=True) + self._status( + f"After clip: {len(annotations)} segments " + f"({n_loaded - len(annotations)} removed outside study area)" + ) + + audible_actors = [] + inaudible_actors = [] + audible_polylines: list[pv.PolyData] = [] + inaudible_polylines: list[pv.PolyData] = [] + n_plotted = 0 + n_skipped = 0 + for _, annot in tqdm( + annotations.iterrows(), + total=len(annotations), + desc="Building annotation lines", + unit="segment", + ): + for line in iter_plot_linestrings(annot["geometry"]): + polyline = self._annotation_polyline(line) + if polyline.n_points < 2: + n_skipped += 1 + continue + n_plotted += 1 + if annot["audible"]: + audible_polylines.append(polyline) + else: + inaudible_polylines.append(polyline) + + audible_actor = self._add_annotation_lines( + audible_polylines, color=self.audible_annotation_color + ) + if audible_actor is not None: + audible_actors.append(audible_actor) + inaudible_actor = self._add_annotation_lines( + inaudible_polylines, color=self.inaudible_annotation_color + ) + if inaudible_actor is not None: + inaudible_actors.append(inaudible_actor) + + if n_plotted == 0: + self._status( + f"No annotation segments plotted ({n_loaded} loaded; " + f"{n_skipped} degenerate lines skipped). " + "Check CRS and study-area overlap." + ) + return + self._status( + f"Plotted {n_plotted} annotation line(s) in " + f"{len(audible_actors) + len(inaudible_actors)} mesh(es) " + f"({len(audible_polylines)} audible, {len(inaudible_polylines)} inaudible; " + f"{n_skipped} skipped)" + ) + + def toggle_audible(flag): + for actor in audible_actors: + actor.SetVisibility(flag) + + def toggle_inaudible(flag): + for actor in inaudible_actors: + actor.SetVisibility(flag) + + self.plotter.add_checkbox_button_widget( + callback=toggle_audible, + value=True, + position=(10, 100), + size=25, + color_on="deepskyblue", + ) + self.plotter.add_checkbox_button_widget( + callback=toggle_inaudible, + value=True, + position=(10, 60), + size=25, + color_on="red", + ) + + def plot_audible_transits(self, audible_transits_pkl: str | None = None) -> None: + if audible_transits_pkl: + print(f"Loading audible transits from {audible_transits_pkl}") + listener = AudibleTransits.from_pickle(audible_transits_pkl) + else: + print("Loading audible transits") + matches = glob.glob( + os.path.join( + self.project_dir, + self.unit + self.site, + "Output_Data", + "AUDIBLE_TRANSITS", + f"3D*{self.year}-01-01*Active Space {self.year}*", + "AudibleTransits_object.pkl", + ) + ) + if len(matches) == 0: + print("No audible transits pkl file found found") + return + listener = AudibleTransits.from_pickle(matches[0]) + + tracks = listener.tracks + if tracks.empty: + print("Audible transits is empty, skipping") + return + + print(f"{len(tracks)} audible transits") + if len(tracks) > self.max_tracks: + print("Too many, sampling") + tracks = tracks.sample(self.max_tracks, random_state=4) + print(f"Showing {len(tracks)} transits") + + tracks = tracks.to_crs(self.crs) + + actors = [] + for _, track in tracks.iterrows(): + polyline = create_polyline_3d(track["interp_geometry"]) + actor = self.plotter.add_mesh( + polyline, color=self.audible_transits_color, point_size=2, line_width=2 + ) + actors.append(actor) + + def toggle(flag): + for actor in actors: + actor.SetVisibility(flag) + + self.plotter.add_checkbox_button_widget( + callback=toggle, + value=True, + position=(10, 180), + size=25, + color_on="purple", + ) + + def plot_vessel_tracks( + self, start_date: str | None = None, end_date: str | None = None + ) -> None: + """Plot MXAK AIS vessel transits at the local sea surface.""" + start_date = start_date or f"{self.year}-01-01" + end_date = end_date or f"{self.year}-12-31" + + try: + ais_path = cfg.read("data", "ais") + except KeyError: + print("No [data] ais path in config, skipping vessel tracks.") + return + + print(f"Querying AIS tracks from {start_date} to {end_date}") + try: + raw_tracks = query_ais_mxak( + ais_path=ais_path, + start_date=start_date, + end_date=end_date, + mask=self.study_area, + ) + except AssertionError as exc: + print(f"No AIS tracks loaded: {exc}") + return + tracks = Tracks(raw_tracks, id_col="event_id", datetime_col="TIME", z_col="altitude") + tracks = tracks.to_crs(self.crs) + + track_ids = tracks["track_id"].drop_duplicates() + print(f"{len(track_ids)} vessel transits ({len(tracks)} points)") + if len(track_ids) > self.max_tracks: + print(f"More than {self.max_tracks}, sampling") + selected = track_ids.sample(self.max_tracks, replace=False, random_state=3) + tracks = tracks[tracks["track_id"].isin(selected)] + print(f"Showing {selected.nunique()} transits") + + actors = [] + for _track_id, group in tracks.groupby("track_id", sort=False): + group = group.sort_values("point_dt") + line = track_points_to_linestring(group) + if line.is_empty or line.length == 0: + continue + polyline = self._sea_surface_polyline(line) + actor = self._add_track_line(polyline, color=self.vessel_track_color) + actors.append(actor) + + if not actors: + print("No vessel tracks to plot.") + return + + def toggle(flag): + for actor in actors: + actor.SetVisibility(flag) + + self.plotter.add_checkbox_button_widget( + callback=toggle, + value=True, + position=(10, 220), + size=25, + color_on=self.vessel_track_color, + ) + + def setup_orientation_widgets(self) -> None: + """Bottom-left E/N/Z axes (+Y = north in UTM).""" + self.plotter.add_axes(**utm_orientation_axes_kwargs()) + + def setup_z_scale(self) -> None: + self.plotter.set_scale(1, 1, 2) + + def toggle_z_scale(flag): + self.plotter.set_scale(1, 1, 2 if flag else 1) + + self.plotter.add_checkbox_button_widget( + callback=toggle_z_scale, + value=True, + position=(10, 140), + size=25, + color_on=self.z_scale_toggle_color, + ) diff --git a/tests/scripts/test_viz.py b/tests/viz/test_viz.py similarity index 95% rename from tests/scripts/test_viz.py rename to tests/viz/test_viz.py index 02cd8db..e88aea3 100644 --- a/tests/scripts/test_viz.py +++ b/tests/viz/test_viz.py @@ -5,7 +5,7 @@ import pytest from shapely.geometry import LineString, MultiLineString, Point -from nps_active_space.scripts.viz import ( +from nps_active_space.viz import ( Visualizer, annotation_z_profile, create_polyline_3d, @@ -115,7 +115,7 @@ def fake_elevation(_dem, _lon, _lat): return 0.0 monkeypatch.setattr( - "nps_active_space.scripts.viz.get_elevation", + "nps_active_space.viz.elevation.get_elevation", fake_elevation, ) dense, z_vals = sea_surface_z_profile( @@ -132,7 +132,7 @@ def test_negative_dem_clamped_to_sea_level(self, monkeypatch): line = LineString([(0, 0), (1, 0)]) monkeypatch.setattr( - "nps_active_space.scripts.viz.get_elevation", + "nps_active_space.viz.elevation.get_elevation", lambda _dem, _lon, _lat: -12.0, ) _, z_vals = sea_surface_z_profile( @@ -263,3 +263,11 @@ def test_builds_line_from_points(self): ) line = track_points_to_linestring(gdf) assert len(line.coords) == 3 + + +class TestScriptsVizShim: + def test_scripts_viz_reexports_visualizer(self): + from nps_active_space.scripts import viz as scripts_viz + + assert scripts_viz.Visualizer is Visualizer + assert scripts_viz.parse_deployment is parse_deployment From ca7ccadbbb85918f6249730b29caa03f75b6aaf7 Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 13:00:10 -0800 Subject: [PATCH 13/17] test --- nps_active_space/viz/assets/waypoint_icon.png | Bin 0 -> 1762 bytes nps_active_space/viz/markers.py | 39 ++++++++++++++++++ nps_active_space/viz/visualizer.py | 9 +++- tests/viz/test_viz.py | 10 +++++ 4 files changed, 56 insertions(+), 2 deletions(-) create mode 100644 nps_active_space/viz/assets/waypoint_icon.png diff --git a/nps_active_space/viz/assets/waypoint_icon.png b/nps_active_space/viz/assets/waypoint_icon.png new file mode 100644 index 0000000000000000000000000000000000000000..805c5b7fae4f3d1eccdd87482d7a3108e93099a3 GIT binary patch literal 1762 zcmV<81|9i{P)V>bE-^4JGb|uzbaZfYIxjD6VRUe8Z***FVlHoT zXD?8q4*#m!Do_4hL7>V*r=(glV;_VHVuB`tCXBKmY%qWt?+#6&TQkU_cjw z0bK|NbRihfgCQ(7697ft!I7bmy(DY(3jmfVvh?F!Cu3GM1;7pf>6EVd=Mp~KdkoomC1|h17N_KdE+tX@nP1ibVZiFb1NCMstx6wA7hNk zTleJSovcF+wIc?>fRC0>MM_*0CF=nUQe>&Yo#)J(M5Od=1GenSbD|l5zULd;e?eo* z6-tHxaA>JHM8>RY4|2}aNTwN4+cdjneY|BPMu#Kep}vqM``tA(U&bdt(eWd&$gFDJ z-F0ZrV~nxtx+WYftkf+a85mZkkH+leVeswVP2|-#8^wUX*ss60zCToqmuHNEcQ>yf z0FwZG>yBlXNROO?a&!j9Jwf9T6l~*Ei{wnGhcZgMo_aBf!TS zNdv=m-&NN)TSd`_t7oR8PP-kguE-sxid`tO{>*6qH!3+iuv>*9OQgDtBnSKHE1W;i z)ZDO*Hmw~6Wn^bwk6`~ERT)k&b(jz8PeQ#VO#Nw=u>d6j`cq2N_aTf6=PgbM$Yl_;@dz z1e<$%o2AIoEVnE(i_h^M$&(VJks2=yM5jOpED_dupckbtT{hv>t^1H$P<~tOkRS8C-U3dIesgu#Npy96jj!tzNtl+H(|m4BAAgJ+YgUKg>|&k&}>An z01*eF@KcJz10pic{7kKF3rwwT+-Pd$8^8LMm!56ll)jbU@zQgA^Ueddd5O-9#S6}Z zA~PA1!7y-8WN%E!_#9&fM?V1|%u&K>z7~Ex6fBuET*!3mJFeKih;+1aZ31od5s;07*qoM6N<$ Eg8Xtavj6}9 literal 0 HcmV?d00001 diff --git a/nps_active_space/viz/markers.py b/nps_active_space/viz/markers.py index 93e18e8..200c029 100644 --- a/nps_active_space/viz/markers.py +++ b/nps_active_space/viz/markers.py @@ -1,5 +1,27 @@ """Scene orientation markers and viewport chrome.""" +from __future__ import annotations + +from pathlib import Path + +import pyvista as pv + +WINDOW_TITLE = "NPS ActiveSpace Visualization" +_ASSETS_DIR = Path(__file__).resolve().parent / "assets" +DEFAULT_WINDOW_ICON = _ASSETS_DIR / "waypoint_icon.png" + + +def default_window_icon_path() -> Path: + """Bundled waypoint/map-pin PNG for the native window icon.""" + return DEFAULT_WINDOW_ICON + + +def apply_window_icon(plotter: pv.Plotter, icon_path: Path | None = None) -> None: + """Apply the window icon when a PNG is available (no-op if missing).""" + path = DEFAULT_WINDOW_ICON if icon_path is None else icon_path + if path.is_file(): + set_window_icon(plotter, path) + def utm_orientation_axes_kwargs() -> dict: """Orientation-marker kwargs for UTM scenes (+X east, +Y north). @@ -15,3 +37,20 @@ def utm_orientation_axes_kwargs() -> dict: "zlabel": "Z", "viewport": (0, 0, 0.2, 0.2), } + + +def set_window_icon(plotter: pv.Plotter, icon_path: str | Path) -> None: + """Set the native window icon from a PNG file. + + Works on Windows and many Linux window managers. macOS often ignores VTK + window icons; use a bundled .icns via a Qt wrapper if dock icon matters. + """ + from vtkmodules.vtkIOImage import vtkPNGReader + + path = Path(icon_path) + if not path.is_file(): + raise FileNotFoundError(path) + reader = vtkPNGReader() + reader.SetFileName(str(path)) + reader.Update() + plotter.render_window.SetIcon(reader.GetOutput()) diff --git a/nps_active_space/viz/visualizer.py b/nps_active_space/viz/visualizer.py index d21ff57..a0f9e0d 100644 --- a/nps_active_space/viz/visualizer.py +++ b/nps_active_space/viz/visualizer.py @@ -26,7 +26,11 @@ ) from nps_active_space.utils.models import Annotations, Tracks from nps_active_space.viz.annotations import format_annotation_summary -from nps_active_space.viz.markers import utm_orientation_axes_kwargs +from nps_active_space.viz.markers import ( + WINDOW_TITLE, + apply_window_icon, + utm_orientation_axes_kwargs, +) from nps_active_space.viz.elevation import ( DemElevationSampler, annotation_z_profile, @@ -90,7 +94,8 @@ def __init__( self.study_area = self.study_area.to_crs(self.crs) self._to_wgs84 = pyproj.Transformer.from_crs(self.crs, "epsg:4326", always_xy=True) - self.plotter = pv.Plotter() + self.plotter = pv.Plotter(title=WINDOW_TITLE) + apply_window_icon(self.plotter) self.plot_dem() self.plot_mic() if do_active: diff --git a/tests/viz/test_viz.py b/tests/viz/test_viz.py index e88aea3..e5c41e4 100644 --- a/tests/viz/test_viz.py +++ b/tests/viz/test_viz.py @@ -197,6 +197,16 @@ def sample_utm_many(self, x, y): class TestOrientationWidgets: + def test_window_title_constant(self): + from nps_active_space.viz.markers import WINDOW_TITLE + + assert WINDOW_TITLE == "NPS ActiveSpace Visualization" + + def test_bundled_waypoint_icon_exists(self): + from nps_active_space.viz.markers import default_window_icon_path + + assert default_window_icon_path().is_file() + def test_utm_axes_use_east_north_up_labels(self): kwargs = utm_orientation_axes_kwargs() assert kwargs["xlabel"] == "E" From 0a48c6249c85311e67301b48fc9abf132d3a54de Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 18 Jun 2026 16:26:20 -0800 Subject: [PATCH 14/17] Some minor path/bugfixes + documentation for --vessels viz option for AIS viz --- .../active_space/layered_active_space.py | 17 ++++-- nps_active_space/scripts/README.md | 7 +++ nps_active_space/scripts/viz.py | 40 -------------- nps_active_space/utils/helpers.py | 18 +++++-- nps_active_space/viz/cli.py | 1 - nps_active_space/viz/visualizer.py | 15 +++++- tests/utils/test_load_layered_activespace.py | 53 +++++++++++++++++++ tests/viz/test_viz.py | 6 +-- 8 files changed, 103 insertions(+), 54 deletions(-) create mode 100644 tests/utils/test_load_layered_activespace.py diff --git a/nps_active_space/active_space/layered_active_space.py b/nps_active_space/active_space/layered_active_space.py index dadd285..c58a30b 100644 --- a/nps_active_space/active_space/layered_active_space.py +++ b/nps_active_space/active_space/layered_active_space.py @@ -22,11 +22,22 @@ def __init__(self, designator, layer_dirs, study_area, gain=None, crs="epsg:4326 self.set_gain(gain) self.fit_pbar = None + if not self.layer_dirs: + raise FileNotFoundError( + f"No active space layer directories provided for {designator}." + ) + # determine min and max gain - import here to avoid circular import from nps_active_space.utils.helpers import omni_to_gain - first_layer_dir = list(self.layer_dirs.values())[0] - active_names = glob.glob(os.path.join(first_layer_dir, "*_O_*.geojson")) - gains = list(map(lambda f: omni_to_gain(f), active_names)) + active_names = [] + for layer_dir in self.layer_dirs.values(): + active_names.extend(glob.glob(os.path.join(layer_dir, "*_O_*.geojson"))) + if not active_names: + raise FileNotFoundError( + f"No *_O_*.geojson active space files found for {designator} in:\n" + + "\n".join(f" {d}" for d in self.layer_dirs.values()) + ) + gains = [omni_to_gain(path) for path in active_names] self.min_gain = min(gains) self.max_gain = max(gains) diff --git a/nps_active_space/scripts/README.md b/nps_active_space/scripts/README.md index 622e9c8..1ff7350 100644 --- a/nps_active_space/scripts/README.md +++ b/nps_active_space/scripts/README.md @@ -458,10 +458,13 @@ This script is used to visualize select geospatial objects relevant to the `nps_ | `-s`, `--active-space` | If included, load and plot the active space | | `-a`, `--annotations` | If included, load and plot annotations | | `-t`, `--audible-transits` | If included, load and plot audible transits | +| `-v`, `--vessels` | If included, load and plot MXAK AIS vessel tracks at sea level | | `--all` | Load and plot all geospatial objects (shorthand for `--active-space --annotations --audible-transits`) | | `-m`, `--max-tracks` | **_default 500_**
Maximum number of annotation tracks or audible transits to show | | `--annotation-file` | **_default to deployment dir_**
Path to .geojson file from which to load annotations | | `--transits-pkl` | **_default to deployment dir_**
Path to .pkl file from which to load audible transits | +| `--start-date` | AIS query start date (YYYY-MM-DD). Default: Jan 1 of deployment year. | +| `--end-date` | AIS query end date (YYYY-MM-DD). Default: Dec 31 of deployment year. | | `--terraced` | If included, render the active space as a terraced surface instead of contours | | `--fill-layers` | If included, fill the interior of each active space contour polygon | @@ -471,6 +474,10 @@ Example executions: $ python -u -W ignore nps_active_space/scripts/viz.py DENATRLA2024 -e production --all ``` +```bash +$ python -u -W ignore nps_active_space/scripts/viz.py GLBALSTL2024 -e production -s -a -v --terraced +``` + ```bash $ python -u -W ignore nps_active_space/scripts/viz.py DENATRLA2024 -e production -g 15.0 -s -a -m 700 --terraced ``` diff --git a/nps_active_space/scripts/viz.py b/nps_active_space/scripts/viz.py index fca8044..cffb903 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -1,46 +1,6 @@ """CLI entry: python -m nps_active_space.scripts.viz or scripts/viz.py path.""" -from nps_active_space.viz import ( - DemElevationSampler, - Visualizer, - annotation_z_profile, - create_polyline_3d, - densify_linestring, - format_annotation_summary, - is_surface_track, - iter_plot_linestrings, - parse_deployment, - parse_existing_file, - parse_iso_date, - parse_max_tracks, - resolve_viz_plot_flags, - sea_surface_z_profile, - track_points_to_linestring, - utm_orientation_axes_kwargs, - vertex_z_from_coord, -) from nps_active_space.viz.cli import main -__all__ = [ - "DemElevationSampler", - "Visualizer", - "annotation_z_profile", - "create_polyline_3d", - "densify_linestring", - "format_annotation_summary", - "is_surface_track", - "iter_plot_linestrings", - "main", - "parse_deployment", - "parse_existing_file", - "parse_iso_date", - "parse_max_tracks", - "resolve_viz_plot_flags", - "sea_surface_z_profile", - "track_points_to_linestring", - "utm_orientation_axes_kwargs", - "vertex_z_from_coord", -] - if __name__ == "__main__": main() diff --git a/nps_active_space/utils/helpers.py b/nps_active_space/utils/helpers.py index f46a8bd..feffaec 100644 --- a/nps_active_space/utils/helpers.py +++ b/nps_active_space/utils/helpers.py @@ -64,14 +64,22 @@ def omni_to_gain(omni_source: str) -> float: def load_layered_activespace(project_dir, unit, site, year, gain=None, crs="epsg:4326"): - prefix = Rf"{project_dir}\{unit}{site}\Output_Data\ACTIVESPACES" + prefix = os.path.join(project_dir, unit + site, "Output_Data", "ACTIVESPACES") layer_dirs = {} - output_dirs = glob.glob(Rf"{prefix}\{unit}{site}{year}_*m") + output_dirs = glob.glob(os.path.join(prefix, f"{unit}{site}{year}_*m")) for dir in output_dirs: + if not glob.glob(os.path.join(dir, "*_O_*.geojson")): + continue altitude = int(os.path.basename(dir).split("_")[1].split("m")[0]) layer_dirs[altitude] = dir + if not layer_dirs: + raise FileNotFoundError( + f"No active space layers with GeoJSON output found for {unit}{site}{year} under:\n" + f" {prefix}\n" + f"Expected directories like {unit}{site}{year}_0m containing *_O_*.geojson files." + ) study_area = load_studyarea(project_dir, unit, site, year) - return LayeredActiveSpace(unit+site+year, layer_dirs, study_area, gain, crs) + return LayeredActiveSpace(f"{unit}{site}{year}", layer_dirs, study_area, gain, crs) def load_activespace(project_dir, unit, site, year, gain, altitude_m=None, crs=None): @@ -106,7 +114,7 @@ def load_activespace(project_dir, unit, site, year, gain, altitude_m=None, crs=N # pick middle altitude if no altitude provided if altitude_m is None: - altitude_dirs = glob.glob(Rf"{prefix}\{unit}{site}{year}_*m") + altitude_dirs = glob.glob(os.path.join(prefix, f"{unit}{site}{year}_*m")) altitudes = [] for dir in altitude_dirs: altitudes.append(int(os.path.basename(dir).split("_")[1].split("m")[0])) @@ -118,7 +126,7 @@ def load_activespace(project_dir, unit, site, year, gain, altitude_m=None, crs=N sign = "-" if gain < 0 else "+" gain_string = str(np.abs(int(10*gain))).zfill(3) usy = f"{unit}{site}{year}" - path = Rf"{prefix}\{usy}_{altitude_m}m\{usy}_O_{sign}{gain_string}.geojson" + path = os.path.join(prefix, f"{usy}_{altitude_m}m", f"{usy}_O_{sign}{gain_string}.geojson") active_space = gpd.read_file(path) if crs is not None: diff --git a/nps_active_space/viz/cli.py b/nps_active_space/viz/cli.py index bc5c773..9b03908 100644 --- a/nps_active_space/viz/cli.py +++ b/nps_active_space/viz/cli.py @@ -162,7 +162,6 @@ def main() -> None: args = parser.parse_args() unit, site, year = args.deployment - print(unit, site, year) do_active, do_annotations, do_transits, do_vessels = resolve_viz_plot_flags( active_space=args.active_space, diff --git a/nps_active_space/viz/visualizer.py b/nps_active_space/viz/visualizer.py index a0f9e0d..d22164a 100644 --- a/nps_active_space/viz/visualizer.py +++ b/nps_active_space/viz/visualizer.py @@ -203,6 +203,14 @@ def plot_activespace(self, terraced: bool = False, gain: float | None = None) -> active_3d = load_layered_activespace( self.project_dir, self.unit, self.site, self.year, gain, self.crs ) + if not active_3d.activespaces: + self._status( + f"No active space geometry loaded for {self.deployment} at {gain}dB; skipping." + ) + return + self._status( + f"Loaded active space layers (m): {sorted(active_3d.activespaces.keys())} at {gain}dB" + ) if terraced: self.plot_terraced_activespace(active_3d) else: @@ -244,7 +252,10 @@ def plot_active_layer(self, active_layer, elevation: float, i: int = 0): for line in active_to_linestrings(active_layer): polyline = create_polyline_3d(line, z=elevation) actor = self.plotter.add_mesh( - polyline, color=self.activespace_color, point_size=2, line_width=2 + polyline, + color=self.activespace_color, + line_width=4, + render_lines_as_tubes=True, ) line_actors.append(actor) @@ -464,7 +475,7 @@ def plot_audible_transits(self, audible_transits_pkl: str | None = None) -> None ) ) if len(matches) == 0: - print("No audible transits pkl file found found") + self._status("No audible transits pkl file found") return listener = AudibleTransits.from_pickle(matches[0]) diff --git a/tests/utils/test_load_layered_activespace.py b/tests/utils/test_load_layered_activespace.py new file mode 100644 index 0000000..5e93b34 --- /dev/null +++ b/tests/utils/test_load_layered_activespace.py @@ -0,0 +1,53 @@ +import pytest +from shapely.geometry import Polygon + +import geopandas as gpd + +from nps_active_space.utils.helpers import load_layered_activespace + + +class TestLoadLayeredActivespace: + def test_finds_layer_dirs_with_posix_paths(self, tmp_path): + project_dir = tmp_path / "project" + site_dir = project_dir / "GLBALSTL" / "Output_Data" / "ACTIVESPACES" + for altitude in (0, 300): + layer_dir = site_dir / f"GLBALSTL2024_{altitude}m" + layer_dir.mkdir(parents=True) + gpd.GeoDataFrame( + geometry=[Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])], + crs="EPSG:4326", + ).to_file(layer_dir / "GLBALSTL2024_O_+100.geojson", driver="GeoJSON") + + study_shp = project_dir / "GLBALSTL" / "GLBALSTL_study_area.shp" + study_shp.parent.mkdir(parents=True, exist_ok=True) + gpd.GeoDataFrame( + geometry=[Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])], + crs="EPSG:4326", + ).to_file(study_shp) + + model = load_layered_activespace(str(project_dir), "GLBA", "LSTL", 2024, gain=10.0) + assert sorted(model.layer_dirs) == [0, 300] + + def test_raises_when_no_layers(self, tmp_path): + project_dir = tmp_path / "project" + (project_dir / "GLBALSTL" / "Output_Data" / "ACTIVESPACES").mkdir(parents=True) + gpd.GeoDataFrame( + geometry=[Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])], + crs="EPSG:4326", + ).to_file(project_dir / "GLBALSTL" / "GLBALSTL_study_area.shp") + + with pytest.raises(FileNotFoundError, match="No active space layers with GeoJSON"): + load_layered_activespace(str(project_dir), "GLBA", "LSTL", 2024) + + def test_ignores_empty_layer_directories(self, tmp_path): + project_dir = tmp_path / "project" + site_dir = project_dir / "GLBALSTL" / "Output_Data" / "ACTIVESPACES" + (site_dir / "GLBALSTL2024_0m").mkdir(parents=True) + (site_dir / "GLBALSTL2024_300m").mkdir(parents=True) + gpd.GeoDataFrame( + geometry=[Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])], + crs="EPSG:4326", + ).to_file(project_dir / "GLBALSTL" / "GLBALSTL_study_area.shp") + + with pytest.raises(FileNotFoundError, match="No active space layers with GeoJSON"): + load_layered_activespace(str(project_dir), "GLBA", "LSTL", 2024) diff --git a/tests/viz/test_viz.py b/tests/viz/test_viz.py index e5c41e4..fa99f73 100644 --- a/tests/viz/test_viz.py +++ b/tests/viz/test_viz.py @@ -276,8 +276,8 @@ def test_builds_line_from_points(self): class TestScriptsVizShim: - def test_scripts_viz_reexports_visualizer(self): + def test_scripts_viz_delegates_to_viz_cli(self): from nps_active_space.scripts import viz as scripts_viz + from nps_active_space.viz.cli import main - assert scripts_viz.Visualizer is Visualizer - assert scripts_viz.parse_deployment is parse_deployment + assert scripts_viz.main is main From e434c284563b49d43f9f43d3cafa2e0750b2ef11 Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Thu, 25 Jun 2026 16:48:44 -0800 Subject: [PATCH 15/17] WIP: viz --track-source and shared load_tracks in utils. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Move track loading to utils for ground truthing and viz; add --track-source CLI with source-specific z handling. Not finalized — needs review and on-network smoke tests. Co-authored-by: Cursor --- .../scripts/run_ground_truthing.py | 2 +- nps_active_space/utils/ais/query.py | 4 +- .../{ground_truthing => utils}/load_tracks.py | 34 ++++---- nps_active_space/viz/cli.py | 26 +++--- nps_active_space/viz/geometry.py | 16 +++- nps_active_space/viz/visualizer.py | 84 ++++++++++++------- tests/viz/test_viz.py | 26 +++++- 7 files changed, 123 insertions(+), 69 deletions(-) rename nps_active_space/{ground_truthing => utils}/load_tracks.py (70%) diff --git a/nps_active_space/scripts/run_ground_truthing.py b/nps_active_space/scripts/run_ground_truthing.py index 872fdee..20fd8da 100644 --- a/nps_active_space/scripts/run_ground_truthing.py +++ b/nps_active_space/scripts/run_ground_truthing.py @@ -2,7 +2,7 @@ import os import iyore import nps_active_space.ground_truthing as app -from nps_active_space.ground_truthing.load_tracks import load_tracks +from nps_active_space.utils.load_tracks import load_tracks from nps_active_space.utils.enums import TrackSource from nps_active_space.utils.models import Nvspl import nps_active_space.utils.config as cfg diff --git a/nps_active_space/utils/ais/query.py b/nps_active_space/utils/ais/query.py index 8b2fd53..d44b5d7 100644 --- a/nps_active_space/utils/ais/query.py +++ b/nps_active_space/utils/ais/query.py @@ -13,7 +13,7 @@ def query_ais_mxak( - ais_path: Path, + ais_path: str | Path, start_date: str, end_date: str, mask: Optional[gpd.GeoDataFrame] = None, @@ -39,6 +39,8 @@ def query_ais_mxak( MxakAis Vessel track points from the MXAK archive. """ + ais_path = Path(ais_path) + if mask is not None: if mask_buffer_distance: ak_albers_mask = mask.to_crs(epsg=3338) diff --git a/nps_active_space/ground_truthing/load_tracks.py b/nps_active_space/utils/load_tracks.py similarity index 70% rename from nps_active_space/ground_truthing/load_tracks.py rename to nps_active_space/utils/load_tracks.py index 14c19db..341253e 100644 --- a/nps_active_space/ground_truthing/load_tracks.py +++ b/nps_active_space/utils/load_tracks.py @@ -4,13 +4,13 @@ import geopandas as gpd import nps_active_space.utils.config as cfg from nps_active_space.utils.ais import query_ais_mxak -from nps_active_space.utils.helpers import create_overflights_engine, query_adsb, query_tracks from nps_active_space.utils.enums import TrackSource +from nps_active_space.utils.helpers import create_overflights_engine, query_adsb, query_tracks from nps_active_space.utils.models import Microphone, Tracks from nps_active_space.utils.time_utils import site_timezone_name, utc_naive_to_site_naive -class GroundTruthingTracks(NamedTuple): +class LoadedTracks(NamedTuple): tracks: Tracks faa_path: str | None faa_corrections_path: str | None @@ -18,8 +18,8 @@ class GroundTruthingTracks(NamedTuple): def _faa_paths() -> tuple[str, str]: return ( - cfg.read('project', 'FAA_Releasable_db'), - cfg.read('project', 'FAA_type_corrections'), + cfg.read("project", "FAA_Releasable_db"), + cfg.read("project", "FAA_type_corrections"), ) @@ -27,33 +27,33 @@ def _load_adsb_tracks( start_date: str, end_date: str, study_area: gpd.GeoDataFrame, -) -> GroundTruthingTracks: +) -> LoadedTracks: raw_tracks = query_adsb( - adsb_path=cfg.read('data', 'adsb'), + adsb_path=cfg.read("data", "adsb"), start_date=start_date, end_date=end_date, mask=study_area, ) - tracks = Tracks(raw_tracks, id_col='flight_id', datetime_col='TIME', z_col='altitude') + tracks = Tracks(raw_tracks, id_col="flight_id", datetime_col="TIME", z_col="altitude") faa_path, faa_corrections_path = _faa_paths() - return GroundTruthingTracks(tracks, faa_path, faa_corrections_path) + return LoadedTracks(tracks, faa_path, faa_corrections_path) def _load_gps_tracks( start_date: str, end_date: str, study_area: gpd.GeoDataFrame, -) -> GroundTruthingTracks: - engine = create_overflights_engine(cfg.read('database:overflights')) +) -> LoadedTracks: + engine = create_overflights_engine(cfg.read("database:overflights")) raw_tracks = query_tracks( engine=engine, start_date=start_date, end_date=end_date, mask=study_area, ) - tracks = Tracks(raw_tracks, 'flight_id', 'ak_datetime', 'altitude_m') + tracks = Tracks(raw_tracks, "flight_id", "ak_datetime", "altitude_m") faa_path, faa_corrections_path = _faa_paths() - return GroundTruthingTracks(tracks, faa_path, faa_corrections_path) + return LoadedTracks(tracks, faa_path, faa_corrections_path) def _load_ais_tracks( @@ -61,17 +61,17 @@ def _load_ais_tracks( end_date: str, study_area: gpd.GeoDataFrame, microphone: Microphone, -) -> GroundTruthingTracks: +) -> LoadedTracks: raw_tracks = query_ais_mxak( ais_path=Path(cfg.read("data", "ais")), start_date=start_date, end_date=end_date, mask=study_area, ) - tracks = Tracks(raw_tracks, id_col='event_id', datetime_col='TIME', z_col='altitude') + tracks = Tracks(raw_tracks, id_col="event_id", datetime_col="TIME", z_col="altitude") site_tz = site_timezone_name(microphone.lat, microphone.lon) tracks["point_dt"] = utc_naive_to_site_naive(tracks["point_dt"], site_tz) - return GroundTruthingTracks(tracks, None, None) + return LoadedTracks(tracks, None, None) def load_tracks( @@ -81,8 +81,8 @@ def load_tracks( end_date: str, study_area: gpd.GeoDataFrame, microphone: Microphone, -) -> GroundTruthingTracks: - """Load flight tracks and FAA lookup paths for a ground-truthing session.""" +) -> LoadedTracks: + """Load tracks for a deployment window and optional FAA lookup paths.""" match source: case TrackSource.ADSB: return _load_adsb_tracks(start_date, end_date, study_area) diff --git a/nps_active_space/viz/cli.py b/nps_active_space/viz/cli.py index 9b03908..387a1b0 100644 --- a/nps_active_space/viz/cli.py +++ b/nps_active_space/viz/cli.py @@ -4,6 +4,7 @@ from datetime import datetime from pathlib import Path +from nps_active_space.utils.enums import TrackSource from nps_active_space.viz.visualizer import Visualizer @@ -65,17 +66,16 @@ def resolve_viz_plot_flags( active_space: bool = False, annotations: bool = False, audible_transits: bool = False, - vessels: bool = False, + track_source: TrackSource | None = None, plot_all: bool = False, annotation_file: str | None = None, transits_pkl: str | None = None, -) -> tuple[bool, bool, bool, bool]: +) -> tuple[bool, bool, bool, TrackSource | None]: """Map CLI flags to Visualizer layer toggles.""" do_active = active_space or plot_all do_annotations = annotations or plot_all or annotation_file is not None do_transits = audible_transits or plot_all or transits_pkl is not None - do_vessels = vessels - return do_active, do_annotations, do_transits, do_vessels + return do_active, do_annotations, do_transits, track_source def main() -> None: @@ -112,10 +112,10 @@ def main() -> None: help="If included, load and plot audible transits", ) parser.add_argument( - "-v", - "--vessels", - action="store_true", - help="If included, load and plot MXAK AIS vessel tracks at sea level", + "--track-source", + type=TrackSource, + choices=list(TrackSource), + help="Load and plot tracks from GPS, ADSB, or AIS.", ) parser.add_argument( "--all", @@ -142,12 +142,12 @@ def main() -> None: parser.add_argument( "--start-date", type=lambda d: parse_iso_date(d, arg_name="--start-date"), - help="AIS query start date (YYYY-MM-DD). Default: Jan 1 of deployment year.", + help="Track query start date (YYYY-MM-DD). Default: Jan 1 of deployment year.", ) parser.add_argument( "--end-date", type=lambda d: parse_iso_date(d, arg_name="--end-date"), - help="AIS query end date (YYYY-MM-DD). Default: Dec 31 of deployment year.", + help="Track query end date (YYYY-MM-DD). Default: Dec 31 of deployment year.", ) parser.add_argument( "--terraced", @@ -163,11 +163,11 @@ def main() -> None: args = parser.parse_args() unit, site, year = args.deployment - do_active, do_annotations, do_transits, do_vessels = resolve_viz_plot_flags( + do_active, do_annotations, do_transits, track_source = resolve_viz_plot_flags( active_space=args.active_space, annotations=args.annotations, audible_transits=args.audible_transits, - vessels=args.vessels, + track_source=args.track_source, plot_all=args.all, annotation_file=args.annotation_file, transits_pkl=args.transits_pkl, @@ -182,7 +182,7 @@ def main() -> None: args.gain, do_annotations, do_transits, - do_vessels, + track_source, args.annotation_file, args.transits_pkl, args.start_date, diff --git a/nps_active_space/viz/geometry.py b/nps_active_space/viz/geometry.py index bd066ff..4ed4777 100644 --- a/nps_active_space/viz/geometry.py +++ b/nps_active_space/viz/geometry.py @@ -96,9 +96,19 @@ def flat_sea_surface_polyline(linestring: LineString, offset_m: float) -> pv.Pol return create_polyline_3d(linestring, z=z) -def track_points_to_linestring(track_points: gpd.GeoDataFrame) -> LineString: - """Connect ordered track points into a 2D line in the plot CRS.""" - coords = [(geom.x, geom.y) for geom in track_points.geometry] +def track_points_to_linestring( + track_points: gpd.GeoDataFrame, + *, + include_z: bool = False, +) -> LineString: + """Connect ordered track points into a line in the plot CRS.""" + if include_z and "z" in track_points.columns: + coords = [ + (geom.x, geom.y, float(z) if np.isfinite(z) else 0.0) + for geom, z in zip(track_points.geometry, track_points["z"], strict=True) + ] + else: + coords = [(geom.x, geom.y) for geom in track_points.geometry] if len(coords) < 2: return LineString(coords) return LineString(coords) diff --git a/nps_active_space/viz/visualizer.py b/nps_active_space/viz/visualizer.py index d22164a..ec5ecd9 100644 --- a/nps_active_space/viz/visualizer.py +++ b/nps_active_space/viz/visualizer.py @@ -14,8 +14,8 @@ import nps_active_space.utils.config as cfg from nps_active_space.scripts.run_audible_transits import AudibleTransits -from nps_active_space.utils.ais import query_ais_mxak from nps_active_space.utils.computation import NMSIM_bbox_utm +from nps_active_space.utils.enums import TrackSource from nps_active_space.utils.helpers import ( get_deployment, get_logger, @@ -24,7 +24,8 @@ load_layered_activespace, load_studyarea, ) -from nps_active_space.utils.models import Annotations, Tracks +from nps_active_space.utils.load_tracks import load_tracks +from nps_active_space.utils.models import Annotations from nps_active_space.viz.annotations import format_annotation_summary from nps_active_space.viz.markers import ( WINDOW_TITLE, @@ -56,6 +57,7 @@ class Visualizer: inaudible_annotation_color = "red" audible_transits_color = "purple" vessel_track_color = "cyan" + flight_track_color = "yellow" z_scale_toggle_color = "black" # button toggling z scale sea_surface_offset_m = 5.0 sea_surface_densify_step_m = 100.0 @@ -70,11 +72,11 @@ def __init__( gain: float | None = None, do_annots: bool = False, do_transits: bool = False, - do_vessels: bool = False, + track_source: TrackSource | None = None, annotation_file: str | None = None, audible_transits_pkl: str | None = None, - vessel_start_date: str | None = None, - vessel_end_date: str | None = None, + track_start_date: str | None = None, + track_end_date: str | None = None, terraced: bool = False, fill_layers: bool = False, max_tracks: int = 1000, @@ -104,8 +106,8 @@ def __init__( self.plot_annotations(annotation_file) if do_transits: self.plot_audible_transits(audible_transits_pkl) - if do_vessels: - self.plot_vessel_tracks(vessel_start_date, vessel_end_date) + if track_source is not None: + self.plot_tracks(track_source, track_start_date, track_end_date) self.plotter.enable_terrain_style() self.setup_z_scale() @@ -512,53 +514,73 @@ def toggle(flag): color_on="purple", ) - def plot_vessel_tracks( - self, start_date: str | None = None, end_date: str | None = None + def _flight_track_polyline(self, linestring: LineString) -> pv.PolyData: + """Build a 3D polyline using stored MSL altitudes from flight tracks.""" + return self._annotation_polyline(linestring) + + def plot_tracks( + self, + source: TrackSource, + start_date: str | None = None, + end_date: str | None = None, ) -> None: - """Plot MXAK AIS vessel transits at the local sea surface.""" + """Plot causal tracks from GPS, ADSB, or MXAK AIS for the study window.""" start_date = start_date or f"{self.year}-01-01" end_date = end_date or f"{self.year}-12-31" + microphone = get_deployment(self.project_dir, self.unit, self.site, self.year) + print(f"Querying {source} tracks from {start_date} to {end_date}") try: - ais_path = cfg.read("data", "ais") - except KeyError: - print("No [data] ais path in config, skipping vessel tracks.") - return - - print(f"Querying AIS tracks from {start_date} to {end_date}") - try: - raw_tracks = query_ais_mxak( - ais_path=ais_path, + loaded = load_tracks( + source, start_date=start_date, end_date=end_date, - mask=self.study_area, + study_area=self.study_area, + microphone=microphone, ) - except AssertionError as exc: - print(f"No AIS tracks loaded: {exc}") + except (KeyError, AssertionError, ValueError) as exc: + print(f"No {source} tracks loaded: {exc}") + return + + tracks = loaded.tracks.to_crs(self.crs) + if tracks.empty: + print(f"No {source} tracks loaded.") return - tracks = Tracks(raw_tracks, id_col="event_id", datetime_col="TIME", z_col="altitude") - tracks = tracks.to_crs(self.crs) track_ids = tracks["track_id"].drop_duplicates() - print(f"{len(track_ids)} vessel transits ({len(tracks)} points)") + print(f"{len(track_ids)} {source} tracks ({len(tracks)} points)") if len(track_ids) > self.max_tracks: print(f"More than {self.max_tracks}, sampling") selected = track_ids.sample(self.max_tracks, replace=False, random_state=3) tracks = tracks[tracks["track_id"].isin(selected)] - print(f"Showing {selected.nunique()} transits") + print(f"Showing {selected.nunique()} tracks") + color = ( + self.vessel_track_color + if source is TrackSource.AIS + else self.flight_track_color + ) actors = [] for _track_id, group in tracks.groupby("track_id", sort=False): group = group.sort_values("point_dt") - line = track_points_to_linestring(group) + line = track_points_to_linestring( + group, + include_z=source is not TrackSource.AIS, + ) if line.is_empty or line.length == 0: continue - polyline = self._sea_surface_polyline(line) - actor = self._add_track_line(polyline, color=self.vessel_track_color) + match source: + case TrackSource.AIS: + polyline = self._sea_surface_polyline(line) + case TrackSource.ADSB | TrackSource.GPS: + polyline = self._flight_track_polyline(line) + case _: + raise ValueError(f"Unknown track source: {source}") + actor = self._add_track_line(polyline, color=color) actors.append(actor) if not actors: - print("No vessel tracks to plot.") + print(f"No {source} tracks to plot.") return def toggle(flag): @@ -570,7 +592,7 @@ def toggle(flag): value=True, position=(10, 220), size=25, - color_on=self.vessel_track_color, + color_on=color, ) def setup_orientation_widgets(self) -> None: diff --git a/tests/viz/test_viz.py b/tests/viz/test_viz.py index fa99f73..029c980 100644 --- a/tests/viz/test_viz.py +++ b/tests/viz/test_viz.py @@ -5,6 +5,7 @@ import pytest from shapely.geometry import LineString, MultiLineString, Point +from nps_active_space.utils.enums import TrackSource from nps_active_space.viz import ( Visualizer, annotation_z_profile, @@ -71,15 +72,23 @@ def test_missing_file(self, tmp_path: Path): class TestResolveVizPlotFlags: def test_annotation_file_implies_annotations(self): flags = resolve_viz_plot_flags(annotation_file="/tmp/a.geojson") - assert flags == (False, True, False, False) + assert flags == (False, True, False, None) def test_transits_pkl_implies_transits(self): flags = resolve_viz_plot_flags(transits_pkl="/tmp/t.pkl") - assert flags == (False, False, True, False) + assert flags == (False, False, True, None) def test_all_enables_standard_layers(self): flags = resolve_viz_plot_flags(plot_all=True) - assert flags == (True, True, True, False) + assert flags == (True, True, True, None) + + def test_track_source_passed_through(self): + flags = resolve_viz_plot_flags(track_source=TrackSource.ADSB) + assert flags == (False, False, False, TrackSource.ADSB) + + def test_ais_track_source(self): + flags = resolve_viz_plot_flags(track_source=TrackSource.AIS) + assert flags == (False, False, False, TrackSource.AIS) class TestParseMaxTracks: @@ -274,6 +283,17 @@ def test_builds_line_from_points(self): line = track_points_to_linestring(gdf) assert len(line.coords) == 3 + def test_include_z_adds_altitude_to_coords(self): + import geopandas as gpd + + gdf = gpd.GeoDataFrame( + {"z": [1000.0, 1500.0, 2000.0]}, + geometry=[Point(0, 0), Point(1, 0), Point(2, 0)], + crs="epsg:32608", + ) + line = track_points_to_linestring(gdf, include_z=True) + assert line.coords[1] == pytest.approx((1.0, 0.0, 1500.0)) + class TestScriptsVizShim: def test_scripts_viz_delegates_to_viz_cli(self): From d963e99d0d1667bbdcac2b010f551485f2bb7ad3 Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Fri, 26 Jun 2026 13:18:34 -0800 Subject: [PATCH 16/17] Add viz --track-source and polish shared load_tracks. Replace deprecated -v/--vessels with --track-source (GPS/ADSB/AIS), route viz through load_tracks with include_faa_paths=False, keep ADSB on the NVSPL clock (no site timezone shift), and add loader tests plus docs. Co-authored-by: Cursor --- example_data/README.md | 21 +- nps_active_space/scripts/README.md | 26 +- .../scripts/downsample_annotations.py | 85 ------ nps_active_space/utils/load_tracks.py | 65 +++- nps_active_space/viz/__init__.py | 2 + nps_active_space/viz/cli.py | 47 ++- nps_active_space/viz/visualizer.py | 37 ++- tests/utils/ais/test_query.py | 2 +- tests/utils/ais/test_reader.py | 2 +- tests/utils/test_load_tracks.py | 286 ++++++++++++++++++ tests/viz/test_viz.py | 246 ++++++++++++++- 11 files changed, 683 insertions(+), 136 deletions(-) delete mode 100644 nps_active_space/scripts/downsample_annotations.py create mode 100644 tests/utils/test_load_tracks.py diff --git a/example_data/README.md b/example_data/README.md index 7d94108..bcdc917 100644 --- a/example_data/README.md +++ b/example_data/README.md @@ -8,7 +8,7 @@ Examples days: Paths in the example configs should parse on Mac/Linux or Windows. However, the paths are relative, so scripts need to be run from the repository root. -**Clock conventions:** MXAK AIS timestamps are UTC-naive; NVSPL and ship-visit times are site-local naive (e.g. GLBALSTL ≈ `America/Juneau`). +**Clock conventions:** MXAK AIS CSV timestamps are UTC-naive; `load_tracks` converts them to site-local naive when a deployment microphone is provided (same clock as NVSPL and ship visits). ADS-B TSV Unix epochs are decoded to naive datetimes and used on the NVSPL clock without a site timezone shift (ground-truthing convention). Legacy pre-2020 ADS-B `.txt` (`EarlyAdsb`) is already naive local. NVSPL and ship-visit times are site-local naive (e.g. GLBALSTL ≈ `America/Juneau`). GPS overflights DB times (`ak_datetime`) are already site-local. ## Layout @@ -51,7 +51,8 @@ Paths in the example configs should parse on Mac/Linux or Windows. However, the | `nvspl_archive/.../NVSPL_DENATRLA_2025_06_23_*.txt` | DENATRLA | Hourly SPL for 2025-06-23 | Offline dev (`DENA_example.config`) | ## Local ground truthing (from repo root) -Running the ground truthing script is a great way to view vessel/flight tracks alongside NVSPL data. The following commands should work out of the box after setting up the project virtual environment. + +Running the ground truthing script is a great way to view vessel/flight tracks alongside NVSPL data. The following commands should work out of the box after setting up the project virtual environment. **GLBALSTL (AIS):** `-e GLBA_example` @@ -68,3 +69,19 @@ python -m nps_active_space.scripts.run_ground_truthing \ ``` To rebuild `faa/MASTER_DENATRLA_20250623_sample.txt` after changing the example ADS-B day, extract rows from the full FAA `MASTER.txt` for the unique `ICAO_address` values in that day's TSV (see `FAAReleasable` in `utils/models.py`). + +## Local viz (from repo root) + +**GLBALSTL (AIS):** `-e GLBA_example` + +```bash +python -m nps_active_space.scripts.viz GLBALSTL2024 -e GLBA_example \ + --track-source AIS --start-date 2024-05-24 --end-date 2024-05-24 -m 100 +``` + +**DENATRLA (ADS-B):** `-e DENA_example` + +```bash +python -m nps_active_space.scripts.viz DENATRLA2025 -e DENA_example \ + --track-source ADSB --start-date 2025-06-23 --end-date 2025-06-23 -m 100 +``` diff --git a/nps_active_space/scripts/README.md b/nps_active_space/scripts/README.md index 1ff7350..ca1f981 100644 --- a/nps_active_space/scripts/README.md +++ b/nps_active_space/scripts/README.md @@ -453,18 +453,18 @@ This script is used to visualize select geospatial objects relevant to the `nps_ | command-line arg | description | | -------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------ | | `deployment` (no flag) | **required.**
The deployment name, e.g., DENACATH2018 | -| `-e`, `--environment` | **required.**
The configuration environment to use. _Ex_: To use `production.config` pass `-e production` | +| `-e`, `--environment` | Config environment name (e.g. `production` for `production.config`). Default: `DENA_streamline` | | `-g`, `--gain` | Active space gain, if not the optimal default found in `fits.csv` | | `-s`, `--active-space` | If included, load and plot the active space | | `-a`, `--annotations` | If included, load and plot annotations | | `-t`, `--audible-transits` | If included, load and plot audible transits | -| `-v`, `--vessels` | If included, load and plot MXAK AIS vessel tracks at sea level | -| `--all` | Load and plot all geospatial objects (shorthand for `--active-space --annotations --audible-transits`) | -| `-m`, `--max-tracks` | **_default 500_**
Maximum number of annotation tracks or audible transits to show | +| `--track-source` | Load and plot tracks from **{GPS, ADSB, AIS}**. Uses the same config paths as ground truthing (`[data]` ais/adsb, overflights DB for GPS). Not included in `--all`. Former `-v`/`--vessels` (AIS only) still works but is deprecated. | +| `--all` | Load and plot active space, annotations, and audible transits (does **not** include `--track-source`) | +| `-m`, `--max-tracks` | **_default 500_**
Maximum number of annotation tracks, audible transits, or causal tracks to show | | `--annotation-file` | **_default to deployment dir_**
Path to .geojson file from which to load annotations | | `--transits-pkl` | **_default to deployment dir_**
Path to .pkl file from which to load audible transits | -| `--start-date` | AIS query start date (YYYY-MM-DD). Default: Jan 1 of deployment year. | -| `--end-date` | AIS query end date (YYYY-MM-DD). Default: Dec 31 of deployment year. | +| `--start-date` | Track query start date (YYYY-MM-DD). **Requires `--track-source`.** Default: Jan 1 of deployment year. | +| `--end-date` | Track query end date (YYYY-MM-DD). **Requires `--track-source`.** Default: Dec 31 of deployment year. | | `--terraced` | If included, render the active space as a terraced surface instead of contours | | `--fill-layers` | If included, fill the interior of each active space contour polygon | @@ -475,13 +475,25 @@ $ python -u -W ignore nps_active_space/scripts/viz.py DENATRLA2024 -e production ``` ```bash -$ python -u -W ignore nps_active_space/scripts/viz.py GLBALSTL2024 -e production -s -a -v --terraced +$ python -u -W ignore nps_active_space/scripts/viz.py GLBALSTL2024 -e production -s -a --track-source AIS --terraced ``` ```bash $ python -u -W ignore nps_active_space/scripts/viz.py DENATRLA2024 -e production -g 15.0 -s -a -m 700 --terraced ``` +```bash +python -m nps_active_space.scripts.viz GLBALSTL2024 -e GLBA_example \ + --track-source AIS --start-date 2024-05-24 --end-date 2024-05-24 -m 100 + +python -m nps_active_space.scripts.viz DENATRLA2025 -e DENA_example \ + --track-source ADSB --start-date 2025-06-23 --end-date 2025-06-23 -m 100 +``` + +**Track plotting vs ground truthing:** viz uses the same `load_tracks` loader but draws raw point sequences (not annotation splines), does not apply clock-drift correction, and defaults to the full deployment year unless `--start-date` / `--end-date` are set. Ground truthing uses the NVSPL archive date span and drift files when present. + +**Note:** In viz, `-t` means audible transits; in ground truthing, `-t` means `--track-source`. + ---- ### Generate Active Space Mesh diff --git a/nps_active_space/scripts/downsample_annotations.py b/nps_active_space/scripts/downsample_annotations.py deleted file mode 100644 index 337ac5d..0000000 --- a/nps_active_space/scripts/downsample_annotations.py +++ /dev/null @@ -1,85 +0,0 @@ -"""Downsample dense annotation GeoJSON (1 Hz splines) for smaller files.""" - -from __future__ import annotations - -import argparse -from pathlib import Path - -import geopandas as gpd -from shapely.geometry import LineString - -from nps_active_space.ground_truthing.segments import ( - ANNOTATION_MAX_VERTICES, - compact_geometry, -) - - -def _vertex_count(geom) -> int: - if isinstance(geom, LineString): - return len(geom.coords) - return 1 - - -def downsample_annotations_file( - input_path: Path, - output_path: Path, - *, - max_vertices: int = ANNOTATION_MAX_VERTICES, -) -> tuple[int, int]: - """Read annotation GeoJSON, compact geometries, write output. Returns (before, after) vertex totals.""" - gdf = gpd.read_file(input_path) - before = sum(_vertex_count(g) for g in gdf.geometry) - gdf = gdf.copy() - gdf["geometry"] = gdf.geometry.map( - lambda g: compact_geometry(g, max_vertices=max_vertices) - ) - after = sum(_vertex_count(g) for g in gdf.geometry) - output_path.parent.mkdir(parents=True, exist_ok=True) - gdf.to_file(output_path, driver="GeoJSON") - return before, after - - -def main() -> None: - parser = argparse.ArgumentParser( - description="Downsample dense ground-truthing annotation GeoJSON." - ) - parser.add_argument("input", type=Path, help="Input .geojson") - parser.add_argument( - "output", - type=Path, - nargs="?", - help="Output .geojson (default: _compact.geojson)", - ) - parser.add_argument( - "--max-vertices", - type=int, - default=ANNOTATION_MAX_VERTICES, - help=f"Max vertices per segment (default: {ANNOTATION_MAX_VERTICES})", - ) - args = parser.parse_args() - - input_path = args.input.resolve() - if not input_path.is_file(): - parser.error(f"input not found: {input_path}") - - output_path = ( - args.output.resolve() - if args.output is not None - else input_path.with_name(f"{input_path.stem}_compact.geojson") - ) - - before_bytes = input_path.stat().st_size - before_verts, after_verts = downsample_annotations_file( - input_path, - output_path, - max_vertices=args.max_vertices, - ) - after_bytes = output_path.stat().st_size - - print(f"Wrote {output_path}") - print(f" vertices: {before_verts:,} -> {after_verts:,}") - print(f" file size: {before_bytes / 1e6:.1f} MB -> {after_bytes / 1e6:.1f} MB") - - -if __name__ == "__main__": - main() diff --git a/nps_active_space/utils/load_tracks.py b/nps_active_space/utils/load_tracks.py index 341253e..986526a 100644 --- a/nps_active_space/utils/load_tracks.py +++ b/nps_active_space/utils/load_tracks.py @@ -6,11 +6,13 @@ from nps_active_space.utils.ais import query_ais_mxak from nps_active_space.utils.enums import TrackSource from nps_active_space.utils.helpers import create_overflights_engine, query_adsb, query_tracks -from nps_active_space.utils.models import Microphone, Tracks +from nps_active_space.utils.models import Adsb, EarlyAdsb, Microphone, Tracks from nps_active_space.utils.time_utils import site_timezone_name, utc_naive_to_site_naive class LoadedTracks(NamedTuple): + """Tracks for a deployment window plus optional FAA lookup paths (aircraft only).""" + tracks: Tracks faa_path: str | None faa_corrections_path: str | None @@ -23,10 +25,20 @@ def _faa_paths() -> tuple[str, str]: ) +def _apply_site_local_timestamps(tracks: Tracks, microphone: Microphone | None) -> None: + """Convert UTC-naive ``point_dt`` to site-local naive using deployment lat/lon.""" + if microphone is None: + return + site_tz = site_timezone_name(microphone.lat, microphone.lon) + tracks["point_dt"] = utc_naive_to_site_naive(tracks["point_dt"], site_tz) + + def _load_adsb_tracks( start_date: str, end_date: str, study_area: gpd.GeoDataFrame, + *, + include_faa_paths: bool, ) -> LoadedTracks: raw_tracks = query_adsb( adsb_path=cfg.read("data", "adsb"), @@ -35,14 +47,22 @@ def _load_adsb_tracks( mask=study_area, ) tracks = Tracks(raw_tracks, id_col="flight_id", datetime_col="TIME", z_col="altitude") - faa_path, faa_corrections_path = _faa_paths() - return LoadedTracks(tracks, faa_path, faa_corrections_path) + # TSV Unix → naive datetime is used on the NVSPL clock as-is (ground-truthing convention). + # Legacy .txt EarlyAdsb is also naive local. Do not apply site timezone shift for ADSB. + if not isinstance(raw_tracks, (Adsb, EarlyAdsb)): + raise TypeError(f"Unexpected ADSB query result type: {type(raw_tracks)!r}") + if include_faa_paths: + faa_path, faa_corrections_path = _faa_paths() + return LoadedTracks(tracks, faa_path, faa_corrections_path) + return LoadedTracks(tracks, None, None) def _load_gps_tracks( start_date: str, end_date: str, study_area: gpd.GeoDataFrame, + *, + include_faa_paths: bool, ) -> LoadedTracks: engine = create_overflights_engine(cfg.read("database:overflights")) raw_tracks = query_tracks( @@ -51,16 +71,18 @@ def _load_gps_tracks( end_date=end_date, mask=study_area, ) - tracks = Tracks(raw_tracks, "flight_id", "ak_datetime", "altitude_m") - faa_path, faa_corrections_path = _faa_paths() - return LoadedTracks(tracks, faa_path, faa_corrections_path) + tracks = Tracks(raw_tracks, id_col="flight_id", datetime_col="ak_datetime", z_col="altitude_m") + if include_faa_paths: + faa_path, faa_corrections_path = _faa_paths() + return LoadedTracks(tracks, faa_path, faa_corrections_path) + return LoadedTracks(tracks, None, None) def _load_ais_tracks( start_date: str, end_date: str, study_area: gpd.GeoDataFrame, - microphone: Microphone, + microphone: Microphone | None, ) -> LoadedTracks: raw_tracks = query_ais_mxak( ais_path=Path(cfg.read("data", "ais")), @@ -69,8 +91,7 @@ def _load_ais_tracks( mask=study_area, ) tracks = Tracks(raw_tracks, id_col="event_id", datetime_col="TIME", z_col="altitude") - site_tz = site_timezone_name(microphone.lat, microphone.lon) - tracks["point_dt"] = utc_naive_to_site_naive(tracks["point_dt"], site_tz) + _apply_site_local_timestamps(tracks, microphone) return LoadedTracks(tracks, None, None) @@ -80,14 +101,32 @@ def load_tracks( start_date: str, end_date: str, study_area: gpd.GeoDataFrame, - microphone: Microphone, + microphone: Microphone | None = None, + include_faa_paths: bool = True, ) -> LoadedTracks: - """Load tracks for a deployment window and optional FAA lookup paths.""" + """Load tracks for a deployment window. + + ADSB TSV and legacy ``.txt`` keep parser naive ``point_dt`` on the NVSPL clock (no site + timezone shift). AIS converts UTC-naive timestamps to site-local naive when + ``microphone`` is provided. GPS ``ak_datetime`` is already site-local from the + overflights database. + + ADSB and GPS return FAA lookup paths when ``include_faa_paths`` is true (the + default for ground truthing). Spatial filtering uses ``study_area`` only; + ``microphone`` is not used for extent. + """ match source: case TrackSource.ADSB: - return _load_adsb_tracks(start_date, end_date, study_area) + return _load_adsb_tracks( + start_date, + end_date, + study_area, + include_faa_paths=include_faa_paths, + ) case TrackSource.GPS: - return _load_gps_tracks(start_date, end_date, study_area) + return _load_gps_tracks( + start_date, end_date, study_area, include_faa_paths=include_faa_paths + ) case TrackSource.AIS: return _load_ais_tracks(start_date, end_date, study_area, microphone) case _: diff --git a/nps_active_space/viz/__init__.py b/nps_active_space/viz/__init__.py index dbf658e..621c6ab 100644 --- a/nps_active_space/viz/__init__.py +++ b/nps_active_space/viz/__init__.py @@ -7,6 +7,7 @@ parse_existing_file, parse_iso_date, parse_max_tracks, + resolve_track_source_args, resolve_viz_plot_flags, ) from nps_active_space.viz.markers import utm_orientation_axes_kwargs @@ -42,6 +43,7 @@ "parse_existing_file", "parse_iso_date", "parse_max_tracks", + "resolve_track_source_args", "resolve_viz_plot_flags", "sea_surface_z_profile", "track_points_to_linestring", diff --git a/nps_active_space/viz/cli.py b/nps_active_space/viz/cli.py index 387a1b0..4226dfd 100644 --- a/nps_active_space/viz/cli.py +++ b/nps_active_space/viz/cli.py @@ -1,6 +1,7 @@ from __future__ import annotations import argparse +import sys from datetime import datetime from pathlib import Path @@ -78,6 +79,33 @@ def resolve_viz_plot_flags( return do_active, do_annotations, do_transits, track_source +def resolve_track_source_args( + args: argparse.Namespace, parser: argparse.ArgumentParser +) -> TrackSource | None: + """Apply deprecated vessel flags and validate track date options.""" + track_source = args.track_source + if args.vessels: + if track_source is not None and track_source is not TrackSource.AIS: + parser.error( + "--vessels is an alias for --track-source AIS; " + "cannot combine with another --track-source" + ) + if track_source is None: + print( + "Warning: -v/--vessels is deprecated; use --track-source AIS", + file=sys.stderr, + ) + track_source = TrackSource.AIS + + if (args.start_date or args.end_date) and track_source is None: + parser.error("--start-date and --end-date require --track-source") + + if args.start_date and args.end_date and args.start_date > args.end_date: + parser.error("--start-date must be on or before --end-date") + + return track_source + + def main() -> None: parser = argparse.ArgumentParser() @@ -115,7 +143,13 @@ def main() -> None: "--track-source", type=TrackSource, choices=list(TrackSource), - help="Load and plot tracks from GPS, ADSB, or AIS.", + help="Load and plot causal tracks (GPS, ADSB, or AIS). Not included in --all.", + ) + parser.add_argument( + "-v", + "--vessels", + action="store_true", + help=argparse.SUPPRESS, ) parser.add_argument( "--all", @@ -127,7 +161,7 @@ def main() -> None: "--max-tracks", type=parse_max_tracks, default=500, - help="Maximum number of annotation tracks or audible transits to show.", + help="Maximum tracks to show (annotations, audible transits, or --track-source).", ) parser.add_argument( "--annotation-file", @@ -142,12 +176,12 @@ def main() -> None: parser.add_argument( "--start-date", type=lambda d: parse_iso_date(d, arg_name="--start-date"), - help="Track query start date (YYYY-MM-DD). Default: Jan 1 of deployment year.", + help="Track query start date (YYYY-MM-DD). Requires --track-source. Default: Jan 1 of deployment year.", ) parser.add_argument( "--end-date", type=lambda d: parse_iso_date(d, arg_name="--end-date"), - help="Track query end date (YYYY-MM-DD). Default: Dec 31 of deployment year.", + help="Track query end date (YYYY-MM-DD). Requires --track-source. Default: Dec 31 of deployment year.", ) parser.add_argument( "--terraced", @@ -163,11 +197,12 @@ def main() -> None: args = parser.parse_args() unit, site, year = args.deployment - do_active, do_annotations, do_transits, track_source = resolve_viz_plot_flags( + track_source = resolve_track_source_args(args, parser) + do_active, do_annotations, do_transits, _ = resolve_viz_plot_flags( active_space=args.active_space, annotations=args.annotations, audible_transits=args.audible_transits, - track_source=args.track_source, + track_source=track_source, plot_all=args.all, annotation_file=args.annotation_file, transits_pkl=args.transits_pkl, diff --git a/nps_active_space/viz/visualizer.py b/nps_active_space/viz/visualizer.py index ec5ecd9..6a99eab 100644 --- a/nps_active_space/viz/visualizer.py +++ b/nps_active_space/viz/visualizer.py @@ -1,6 +1,7 @@ from __future__ import annotations import glob +import configparser import os import geopandas as gpd @@ -514,22 +515,23 @@ def toggle(flag): color_on="purple", ) - def _flight_track_polyline(self, linestring: LineString) -> pv.PolyData: - """Build a 3D polyline using stored MSL altitudes from flight tracks.""" - return self._annotation_polyline(linestring) - def plot_tracks( self, source: TrackSource, start_date: str | None = None, end_date: str | None = None, ) -> None: - """Plot causal tracks from GPS, ADSB, or MXAK AIS for the study window.""" + """Plot causal tracks from GPS, ADSB, or MXAK AIS for the study window. + + Uses the same ``load_tracks`` loader as ground truthing but draws raw point + sequences (not annotation splines) and does not apply clock-drift correction. + Prefer explicit ``--start-date`` / ``--end-date``; default is the full year. + """ start_date = start_date or f"{self.year}-01-01" end_date = end_date or f"{self.year}-12-31" microphone = get_deployment(self.project_dir, self.unit, self.site, self.year) - print(f"Querying {source} tracks from {start_date} to {end_date}") + self._status(f"Querying {source} tracks from {start_date} to {end_date}") try: loaded = load_tracks( source, @@ -537,23 +539,30 @@ def plot_tracks( end_date=end_date, study_area=self.study_area, microphone=microphone, + include_faa_paths=False, ) - except (KeyError, AssertionError, ValueError) as exc: - print(f"No {source} tracks loaded: {exc}") + except (configparser.NoSectionError, configparser.NoOptionError) as exc: + self._status(f"No {source} tracks loaded: missing config option {exc!r}") + return + except KeyError as exc: + self._status(f"No {source} tracks loaded: {exc}") + return + except (AssertionError, ValueError) as exc: + self._status(f"No {source} tracks loaded: {exc}") return tracks = loaded.tracks.to_crs(self.crs) if tracks.empty: - print(f"No {source} tracks loaded.") + self._status(f"No {source} tracks loaded.") return track_ids = tracks["track_id"].drop_duplicates() - print(f"{len(track_ids)} {source} tracks ({len(tracks)} points)") + self._status(f"{len(track_ids)} {source} tracks ({len(tracks)} points)") if len(track_ids) > self.max_tracks: - print(f"More than {self.max_tracks}, sampling") + self._status(f"More than {self.max_tracks}, sampling") selected = track_ids.sample(self.max_tracks, replace=False, random_state=3) tracks = tracks[tracks["track_id"].isin(selected)] - print(f"Showing {selected.nunique()} tracks") + self._status(f"Showing {selected.nunique()} tracks") color = ( self.vessel_track_color @@ -573,14 +582,14 @@ def plot_tracks( case TrackSource.AIS: polyline = self._sea_surface_polyline(line) case TrackSource.ADSB | TrackSource.GPS: - polyline = self._flight_track_polyline(line) + polyline = self._annotation_polyline(line) case _: raise ValueError(f"Unknown track source: {source}") actor = self._add_track_line(polyline, color=color) actors.append(actor) if not actors: - print(f"No {source} tracks to plot.") + self._status(f"No {source} tracks to plot.") return def toggle(flag): diff --git a/tests/utils/ais/test_query.py b/tests/utils/ais/test_query.py index 26895bd..34156c1 100644 --- a/tests/utils/ais/test_query.py +++ b/tests/utils/ais/test_query.py @@ -9,7 +9,7 @@ from nps_active_space.utils.ais.reader import MxakAis REPO = Path(__file__).resolve().parents[3] -AIS_FIXTURE = REPO / "example_data" / "AIS" / "MXAK-AIS-GLBA-20250107.csv" +AIS_FIXTURE = REPO / "example_data" / "MXAK-AIS-GLBA" / "MXAK-AIS-GLBA-20250107.csv" AIS_DIR = AIS_FIXTURE.parent METADATA_COLS = ["TIME", "MMSI", "lat", "lon", "ship_name", "shiptype", "altitude", "event_id"] diff --git a/tests/utils/ais/test_reader.py b/tests/utils/ais/test_reader.py index 5b6a71d..f8c02fc 100644 --- a/tests/utils/ais/test_reader.py +++ b/tests/utils/ais/test_reader.py @@ -8,7 +8,7 @@ from nps_active_space.utils.ais.reader import MxakAis REPO = Path(__file__).resolve().parents[3] -AIS_FIXTURE = REPO / "example_data" / "AIS" / "MXAK-AIS-GLBA-20250107.csv" +AIS_FIXTURE = REPO / "example_data" / "MXAK-AIS-GLBA" / "MXAK-AIS-GLBA-20250107.csv" ANNA_MMSI = 368018710 METADATA_COLS = ["TIME", "MMSI", "lat", "lon", "ship_name", "shiptype", "altitude", "event_id"] diff --git a/tests/utils/test_load_tracks.py b/tests/utils/test_load_tracks.py new file mode 100644 index 0000000..64bb6ef --- /dev/null +++ b/tests/utils/test_load_tracks.py @@ -0,0 +1,286 @@ +from __future__ import annotations + +from unittest.mock import MagicMock + +import geopandas as gpd +import pandas as pd +import pytest +from pandas.testing import assert_series_equal +from shapely.geometry import box + +import nps_active_space.utils.config as cfg +import nps_active_space.utils.load_tracks as load_tracks_mod +from nps_active_space.utils.enums import TrackSource +from nps_active_space.utils.load_tracks import LoadedTracks, load_tracks +from nps_active_space.utils.models import Adsb, EarlyAdsb, Microphone, Tracks +from nps_active_space.utils.time_utils import site_timezone_name, utc_naive_to_site_naive + +UTM8 = "epsg:32608" +START = "2024-05-24" +END = "2024-05-24" + + +@pytest.fixture +def study_area() -> gpd.GeoDataFrame: + return gpd.GeoDataFrame(geometry=[box(500_000, 6_000_000, 510_000, 6_010_000)], crs=UTM8) + + +@pytest.fixture +def microphone() -> Microphone: + return Microphone("GLBALSTL2024", lat=58.4, lon=-136.0, z=12.0) + + +def _points_gdf( + *, + id_col: str, + id_value: str, + time_col: str, + z_col: str, + z_value: float, +) -> gpd.GeoDataFrame: + return gpd.GeoDataFrame( + { + id_col: [id_value, id_value], + time_col: [ + pd.Timestamp("2024-05-24 10:00:00"), + pd.Timestamp("2024-05-24 10:05:00"), + ], + z_col: [z_value, z_value + 100.0], + }, + geometry=gpd.points_from_xy([500_100.0, 500_200.0], [6_000_100.0, 6_000_200.0]), + crs=UTM8, + ) + + +def _faa_config_read(section: str, option: str | None = None) -> str: + if section == "project" and option == "FAA_Releasable_db": + return "/faa/MASTER.txt" + if section == "project" and option == "FAA_type_corrections": + return "/faa/corrections.json" + raise KeyError(f"{section}.{option}") + + +class TestLoadTracks: + def test_adsb_passes_query_kwargs_and_normalizes_columns( + self, monkeypatch: pytest.MonkeyPatch, study_area: gpd.GeoDataFrame, microphone: Microphone + ) -> None: + raw = _points_gdf( + id_col="flight_id", + id_value="abc123", + time_col="TIME", + z_col="altitude", + z_value=5000.0, + ) + calls: list[dict] = [] + + def fake_read(section: str, option: str | None = None) -> str: + if section == "data" and option == "adsb": + return "/data/adsb" + return _faa_config_read(section, option) + + def fake_query_adsb(**kwargs) -> Adsb: + calls.append(kwargs) + adsb = Adsb.__new__(Adsb) + gpd.GeoDataFrame.__init__(adsb, raw) + return adsb + + monkeypatch.setattr(cfg, "read", fake_read) + monkeypatch.setattr(load_tracks_mod, "query_adsb", fake_query_adsb) + + loaded = load_tracks( + TrackSource.ADSB, + start_date=START, + end_date=END, + study_area=study_area, + microphone=microphone, + ) + + assert calls == [ + { + "adsb_path": "/data/adsb", + "start_date": START, + "end_date": END, + "mask": study_area, + } + ] + assert isinstance(loaded.tracks, Tracks) + assert loaded.faa_path == "/faa/MASTER.txt" + assert loaded.tracks["track_id"].tolist() == ["abc123", "abc123"] + assert list(loaded.tracks["z"]) == pytest.approx([5000.0, 5100.0]) + assert_series_equal(loaded.tracks["point_dt"], raw["TIME"], check_names=False) + + def test_adsb_early_format_keeps_parser_timestamps( + self, monkeypatch: pytest.MonkeyPatch, study_area: gpd.GeoDataFrame, microphone: Microphone + ) -> None: + raw = _points_gdf( + id_col="flight_id", + id_value="legacy1", + time_col="TIME", + z_col="altitude", + z_value=5000.0, + ) + local_times = raw["TIME"].copy() + + def fake_read(section: str, option: str | None = None) -> str: + if section == "data" and option == "adsb": + return "/data/adsb" + return _faa_config_read(section, option) + + def fake_query_adsb(**kwargs) -> EarlyAdsb: + early = EarlyAdsb.__new__(EarlyAdsb) + gpd.GeoDataFrame.__init__(early, raw) + return early + + monkeypatch.setattr(cfg, "read", fake_read) + monkeypatch.setattr(load_tracks_mod, "query_adsb", fake_query_adsb) + + loaded = load_tracks( + TrackSource.ADSB, + start_date=START, + end_date=END, + study_area=study_area, + microphone=microphone, + ) + + assert_series_equal(loaded.tracks["point_dt"], local_times, check_names=False) + + def test_gps_passes_query_kwargs( + self, monkeypatch: pytest.MonkeyPatch, study_area: gpd.GeoDataFrame, microphone: Microphone + ) -> None: + raw = _points_gdf( + id_col="flight_id", + id_value="gps1", + time_col="ak_datetime", + z_col="altitude_m", + z_value=8000.0, + ) + engine = MagicMock() + track_calls: list[dict] = [] + + def fake_read(section: str, option: str | None = None): + if section == "database:overflights" and option is None: + return { + "name": "db", + "host": "localhost", + "username": "u", + "password": "p", + "port": 5432, + } + return _faa_config_read(section, option) + + def fake_query_tracks(**kwargs) -> gpd.GeoDataFrame: + track_calls.append(kwargs) + return raw + + monkeypatch.setattr(cfg, "read", fake_read) + monkeypatch.setattr(load_tracks_mod, "create_overflights_engine", lambda _cfg: engine) + monkeypatch.setattr(load_tracks_mod, "query_tracks", fake_query_tracks) + + loaded = load_tracks( + TrackSource.GPS, + start_date=START, + end_date=END, + study_area=study_area, + microphone=microphone, + ) + + assert track_calls[0]["engine"] is engine + assert track_calls[0]["start_date"] == START + assert track_calls[0]["mask"] is study_area + assert loaded.tracks["track_id"].iloc[0] == "gps1" + assert loaded.faa_path == "/faa/MASTER.txt" + assert loaded.faa_corrections_path == "/faa/corrections.json" + + def test_adsb_skips_faa_paths_when_disabled( + self, monkeypatch: pytest.MonkeyPatch, study_area: gpd.GeoDataFrame + ) -> None: + raw = _points_gdf( + id_col="flight_id", + id_value="abc123", + time_col="TIME", + z_col="altitude", + z_value=5000.0, + ) + + def fake_read(section: str, option: str | None = None) -> str: + if section == "data" and option == "adsb": + return "/data/adsb" + raise KeyError(f"{section}.{option}") + + def fake_query_adsb(**kwargs) -> Adsb: + adsb = Adsb.__new__(Adsb) + gpd.GeoDataFrame.__init__(adsb, raw) + return adsb + + monkeypatch.setattr(cfg, "read", fake_read) + monkeypatch.setattr(load_tracks_mod, "query_adsb", fake_query_adsb) + + loaded = load_tracks( + TrackSource.ADSB, + start_date=START, + end_date=END, + study_area=study_area, + include_faa_paths=False, + ) + + assert loaded.faa_path is None + assert loaded.faa_corrections_path is None + + def test_ais_without_microphone_keeps_utc_timestamps( + self, monkeypatch: pytest.MonkeyPatch, study_area: gpd.GeoDataFrame + ) -> None: + raw = _points_gdf( + id_col="event_id", + id_value="evt1", + time_col="TIME", + z_col="altitude", + z_value=0.0, + ) + utc_times = raw["TIME"].copy() + + monkeypatch.setattr(cfg, "read", lambda section, option=None: "/data/ais") + monkeypatch.setattr(load_tracks_mod, "query_ais_mxak", lambda **kwargs: raw) + + loaded = load_tracks( + TrackSource.AIS, + start_date=START, + end_date=END, + study_area=study_area, + ) + + assert_series_equal(loaded.tracks["point_dt"], utc_times, check_names=False) + + def test_ais_converts_timestamps_to_site_local( + self, monkeypatch: pytest.MonkeyPatch, study_area: gpd.GeoDataFrame, microphone: Microphone + ) -> None: + raw = _points_gdf( + id_col="event_id", + id_value="evt1", + time_col="TIME", + z_col="altitude", + z_value=0.0, + ) + utc_times = raw["TIME"].copy() + ais_calls: list[dict] = [] + + monkeypatch.setattr(cfg, "read", lambda section, option=None: "/data/ais") + monkeypatch.setattr( + load_tracks_mod, + "query_ais_mxak", + lambda **kwargs: ais_calls.append(kwargs) or raw, + ) + + loaded = load_tracks( + TrackSource.AIS, + start_date=START, + end_date=END, + study_area=study_area, + microphone=microphone, + ) + + assert ais_calls[0]["ais_path"].as_posix() == "/data/ais" + assert ais_calls[0]["mask"] is study_area + site_tz = site_timezone_name(microphone.lat, microphone.lon) + expected = utc_naive_to_site_naive(utc_times, site_tz) + assert_series_equal(loaded.tracks["point_dt"], expected, check_names=False) + assert loaded.faa_path is None diff --git a/tests/viz/test_viz.py b/tests/viz/test_viz.py index 029c980..6a97bde 100644 --- a/tests/viz/test_viz.py +++ b/tests/viz/test_viz.py @@ -1,11 +1,17 @@ import argparse +from configparser import NoOptionError from pathlib import Path +from unittest.mock import MagicMock +import geopandas as gpd import numpy as np +import pandas as pd import pytest -from shapely.geometry import LineString, MultiLineString, Point +from shapely.geometry import LineString, MultiLineString, Point, box from nps_active_space.utils.enums import TrackSource +from nps_active_space.utils.load_tracks import LoadedTracks +from nps_active_space.utils.models import Tracks from nps_active_space.viz import ( Visualizer, annotation_z_profile, @@ -18,6 +24,7 @@ parse_existing_file, parse_iso_date, parse_max_tracks, + resolve_track_source_args, resolve_viz_plot_flags, sea_surface_z_profile, track_points_to_linestring, @@ -82,13 +89,76 @@ def test_all_enables_standard_layers(self): flags = resolve_viz_plot_flags(plot_all=True) assert flags == (True, True, True, None) - def test_track_source_passed_through(self): - flags = resolve_viz_plot_flags(track_source=TrackSource.ADSB) - assert flags == (False, False, False, TrackSource.ADSB) + @pytest.mark.parametrize("source", list(TrackSource)) + def test_track_source_passed_through(self, source: TrackSource) -> None: + flags = resolve_viz_plot_flags(track_source=source) + assert flags == (False, False, False, source) - def test_ais_track_source(self): - flags = resolve_viz_plot_flags(track_source=TrackSource.AIS) - assert flags == (False, False, False, TrackSource.AIS) + +class TestResolveTrackSourceArgs: + def test_dates_require_track_source(self) -> None: + parser = argparse.ArgumentParser() + args = argparse.Namespace( + track_source=None, + vessels=False, + start_date="2024-05-24", + end_date=None, + ) + with pytest.raises(SystemExit): + resolve_track_source_args(args, parser) + + def test_vessels_sets_ais_track_source(self, capsys: pytest.CaptureFixture[str]) -> None: + parser = argparse.ArgumentParser() + args = argparse.Namespace( + track_source=None, + vessels=True, + start_date=None, + end_date=None, + ) + assert resolve_track_source_args(args, parser) is TrackSource.AIS + assert "deprecated" in capsys.readouterr().err + + def test_rejects_vessels_with_other_track_source(self) -> None: + parser = argparse.ArgumentParser() + args = argparse.Namespace( + track_source=TrackSource.ADSB, + vessels=True, + start_date=None, + end_date=None, + ) + with pytest.raises(SystemExit): + resolve_track_source_args(args, parser) + + def test_rejects_inverted_date_range(self) -> None: + parser = argparse.ArgumentParser() + args = argparse.Namespace( + track_source=TrackSource.AIS, + vessels=False, + start_date="2024-05-25", + end_date="2024-05-24", + ) + with pytest.raises(SystemExit): + resolve_track_source_args(args, parser) + + def test_passes_through_track_source_with_dates(self) -> None: + parser = argparse.ArgumentParser() + args = argparse.Namespace( + track_source=TrackSource.ADSB, + vessels=False, + start_date="2024-05-24", + end_date="2024-05-25", + ) + assert resolve_track_source_args(args, parser) is TrackSource.ADSB + + def test_allows_vessels_with_ais_track_source(self) -> None: + parser = argparse.ArgumentParser() + args = argparse.Namespace( + track_source=TrackSource.AIS, + vessels=True, + start_date=None, + end_date=None, + ) + assert resolve_track_source_args(args, parser) is TrackSource.AIS class TestParseMaxTracks: @@ -295,6 +365,168 @@ def test_include_z_adds_altitude_to_coords(self): assert line.coords[1] == pytest.approx((1.0, 0.0, 1500.0)) +class TestPlotTracksSourceRouting: + @staticmethod + def _stub_visualizer() -> Visualizer: + vis = Visualizer.__new__(Visualizer) + vis.unit, vis.site, vis.year = "GLBA", "LSTL", 2024 + vis.project_dir = "/proj" + vis.max_tracks = 500 + vis.crs = "epsg:32608" + vis.study_area = gpd.GeoDataFrame( + geometry=[box(500_000, 6_000_000, 510_000, 6_010_000)], + crs=vis.crs, + ) + vis.plotter = MagicMock() + vis.logger = MagicMock() + vis._dem_sampler = MagicMock() + vis.dem = MagicMock() + vis.sea_surface_offset_m = 5.0 + vis.sea_surface_densify_step_m = 100.0 + vis.vessel_track_color = "cyan" + vis.flight_track_color = "yellow" + return vis + + @staticmethod + def _loaded_tracks(*, flight: bool) -> LoadedTracks: + id_col = "flight_id" if flight else "event_id" + z = 5000.0 if flight else 0.0 + gdf = gpd.GeoDataFrame( + { + id_col: ["t1", "t1"], + "TIME": [ + pd.Timestamp("2024-05-24 10:00"), + pd.Timestamp("2024-05-24 10:05"), + ], + "altitude": [z, z], + }, + geometry=gpd.points_from_xy([500_100.0, 500_200.0], [6_000_100.0, 6_000_200.0]), + crs="epsg:32608", + ) + tracks = Tracks(gdf, id_col=id_col, datetime_col="TIME", z_col="altitude") + return LoadedTracks(tracks, None, None) + + @pytest.mark.parametrize( + ("source", "expect_sea", "expect_flight", "use_flight_tracks"), + [ + (TrackSource.AIS, True, False, False), + (TrackSource.ADSB, False, True, True), + (TrackSource.GPS, False, True, True), + ], + ) + def test_routes_polyline_builder_by_source( + self, + monkeypatch: pytest.MonkeyPatch, + source: TrackSource, + expect_sea: bool, + expect_flight: bool, + use_flight_tracks: bool, + ) -> None: + vis = self._stub_visualizer() + captured_lines: list[LineString] = [] + sea = MagicMock(side_effect=lambda line: captured_lines.append(line) or MagicMock()) + flight = MagicMock(side_effect=lambda line: captured_lines.append(line) or MagicMock()) + monkeypatch.setattr(vis, "_sea_surface_polyline", sea) + monkeypatch.setattr(vis, "_annotation_polyline", flight) + monkeypatch.setattr(vis, "_add_track_line", lambda polyline, *, color: MagicMock()) + monkeypatch.setattr( + "nps_active_space.viz.visualizer.get_deployment", + lambda *args, **kwargs: MagicMock(lat=58.4, lon=-136.0), + ) + monkeypatch.setattr( + "nps_active_space.viz.visualizer.load_tracks", + lambda *args, **kwargs: self._loaded_tracks(flight=use_flight_tracks), + ) + + vis.plot_tracks(source, "2024-05-24", "2024-05-24") + + assert sea.called is expect_sea + assert flight.called is expect_flight + coord_dim = len(captured_lines[0].coords[0]) + if use_flight_tracks: + assert coord_dim == 3 + else: + assert coord_dim == 2 + + def test_forwards_load_tracks_kwargs(self, monkeypatch: pytest.MonkeyPatch) -> None: + vis = self._stub_visualizer() + calls: list[dict] = [] + monkeypatch.setattr(vis, "_add_track_line", lambda polyline, *, color: MagicMock()) + monkeypatch.setattr( + "nps_active_space.viz.visualizer.get_deployment", + lambda *args, **kwargs: MagicMock(lat=58.4, lon=-136.0), + ) + monkeypatch.setattr( + "nps_active_space.viz.visualizer.load_tracks", + lambda *args, **kwargs: calls.append(kwargs) or self._loaded_tracks(flight=True), + ) + + vis.plot_tracks(TrackSource.ADSB, "2024-05-24", "2024-05-25") + + assert calls[0]["start_date"] == "2024-05-24" + assert calls[0]["end_date"] == "2024-05-25" + assert calls[0]["include_faa_paths"] is False + assert calls[0]["study_area"] is vis.study_area + + def test_defaults_to_deployment_year_dates(self, monkeypatch: pytest.MonkeyPatch) -> None: + vis = self._stub_visualizer() + calls: list[dict] = [] + monkeypatch.setattr(vis, "_add_track_line", lambda polyline, *, color: MagicMock()) + monkeypatch.setattr( + "nps_active_space.viz.visualizer.get_deployment", + lambda *args, **kwargs: MagicMock(lat=58.4, lon=-136.0), + ) + monkeypatch.setattr( + "nps_active_space.viz.visualizer.load_tracks", + lambda *args, **kwargs: calls.append(kwargs) or self._loaded_tracks(flight=True), + ) + + vis.plot_tracks(TrackSource.ADSB, None, None) + + assert calls[0]["start_date"] == "2024-01-01" + assert calls[0]["end_date"] == "2024-12-31" + + def test_missing_config_reports_config_error(self, monkeypatch: pytest.MonkeyPatch) -> None: + vis = self._stub_visualizer() + messages: list[str] = [] + vis._status = lambda msg: messages.append(msg) + monkeypatch.setattr( + "nps_active_space.viz.visualizer.get_deployment", + lambda *args, **kwargs: MagicMock(lat=58.4, lon=-136.0), + ) + monkeypatch.setattr( + "nps_active_space.viz.visualizer.load_tracks", + lambda *args, **kwargs: (_ for _ in ()).throw(NoOptionError("adsb", "data")), + ) + + vis.plot_tracks(TrackSource.ADSB, "2024-05-24", "2024-05-24") + + assert any("missing config option" in msg for msg in messages) + + def test_parse_key_error_is_not_labeled_missing_config( + self, monkeypatch: pytest.MonkeyPatch + ) -> None: + vis = self._stub_visualizer() + messages: list[str] = [] + vis._status = lambda msg: messages.append(msg) + monkeypatch.setattr( + "nps_active_space.viz.visualizer.get_deployment", + lambda *args, **kwargs: MagicMock(lat=58.4, lon=-136.0), + ) + monkeypatch.setattr( + "nps_active_space.viz.visualizer.load_tracks", + lambda *args, **kwargs: (_ for _ in ()).throw( + KeyError("MXAK AIS CSV lacks 'TIME' column") + ), + ) + + vis.plot_tracks(TrackSource.AIS, "2024-05-24", "2024-05-24") + + assert messages[-1].startswith("No AIS tracks loaded:") + assert "TIME" in messages[-1] + assert not any("missing config option" in msg for msg in messages) + + class TestScriptsVizShim: def test_scripts_viz_delegates_to_viz_cli(self): from nps_active_space.scripts import viz as scripts_viz From c68e03a0a3155e814db709c4796110f081f7b1d8 Mon Sep 17 00:00:00 2001 From: Elliott Ruebush Date: Mon, 6 Jul 2026 18:52:24 -0400 Subject: [PATCH 17/17] retrigger CI