diff --git a/example_data/README.md b/example_data/README.md index 7d941082..bcdc9174 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/active_space/layered_active_space.py b/nps_active_space/active_space/layered_active_space.py index dadd2853..c58a30b9 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/ground_truthing/load_tracks.py b/nps_active_space/ground_truthing/load_tracks.py deleted file mode 100644 index 14c19db6..00000000 --- a/nps_active_space/ground_truthing/load_tracks.py +++ /dev/null @@ -1,94 +0,0 @@ -from pathlib import Path -from typing import NamedTuple - -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.models import Microphone, Tracks -from nps_active_space.utils.time_utils import site_timezone_name, utc_naive_to_site_naive - - -class GroundTruthingTracks(NamedTuple): - tracks: Tracks - faa_path: str | None - faa_corrections_path: str | None - - -def _faa_paths() -> tuple[str, str]: - return ( - cfg.read('project', 'FAA_Releasable_db'), - cfg.read('project', 'FAA_type_corrections'), - ) - - -def _load_adsb_tracks( - start_date: str, - end_date: str, - study_area: gpd.GeoDataFrame, -) -> GroundTruthingTracks: - raw_tracks = query_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') - faa_path, faa_corrections_path = _faa_paths() - return GroundTruthingTracks(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')) - 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') - faa_path, faa_corrections_path = _faa_paths() - return GroundTruthingTracks(tracks, faa_path, faa_corrections_path) - - -def _load_ais_tracks( - start_date: str, - end_date: str, - study_area: gpd.GeoDataFrame, - microphone: Microphone, -) -> GroundTruthingTracks: - 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') - 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) - - -def load_tracks( - source: TrackSource, - *, - start_date: str, - end_date: str, - study_area: gpd.GeoDataFrame, - microphone: Microphone, -) -> GroundTruthingTracks: - """Load flight tracks and FAA lookup paths for a ground-truthing session.""" - match source: - case TrackSource.ADSB: - return _load_adsb_tracks(start_date, end_date, study_area) - case TrackSource.GPS: - return _load_gps_tracks(start_date, end_date, study_area) - case TrackSource.AIS: - return _load_ais_tracks(start_date, end_date, study_area, microphone) - case _: - raise ValueError(f"Unknown track source: {source}") diff --git a/nps_active_space/ground_truthing/segments.py b/nps_active_space/ground_truthing/segments.py index 3c36a084..0a81a5f5 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/README.md b/nps_active_space/scripts/README.md index 622e9c8e..ca1f981b 100644 --- a/nps_active_space/scripts/README.md +++ b/nps_active_space/scripts/README.md @@ -453,15 +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 | -| `--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` | 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 | @@ -471,10 +474,26 @@ 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 --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/run_ground_truthing.py b/nps_active_space/scripts/run_ground_truthing.py index 872fdee0..20fd8da2 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/scripts/viz.py b/nps_active_space/scripts/viz.py index 6abc109c..cffb9032 100644 --- a/nps_active_space/scripts/viz.py +++ b/nps_active_space/scripts/viz.py @@ -1,438 +1,6 @@ -import os -import numpy as np -import pandas as pd -import geopandas as gpd -import glob -import pyvista as pv -import rasterio -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.scripts.run_audible_transits import AudibleTransits -from nps_active_space.utils.models import Annotations -from nps_active_space.utils.computation import NMSIM_bbox_utm -import nps_active_space.utils.config as cfg -import argparse - -# helper functions ============================================ - -def polygon_to_mesh(polygon, z): - 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): - 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): - geometry = active.boundary.iloc[0] - if isinstance(geometry, MultiLineString): - linestrings = geometry.geoms - elif isinstance(geometry, LineString): - linestrings = [geometry] - return linestrings - - -def create_polyline_3d(linestring, z=None): - coords = np.array(linestring.coords) - if z is not None: - coords = np.column_stack((coords, np.full(coords.shape[0], z))) - assert coords.shape[1] == 3 - 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) - poly.lines = lines - return poly - - - - -# main class ===================================================== - -class Visualizer(): - # color config - activespace_color = "orange" - mic_color = "white" - audible_annotation_color = "deepskyblue" - inaudible_annotation_color = "red" - audible_transits_color = "purple" - z_scale_toggle_color = "black" # button toggling z scale - - def __init__(self, unit, site, year, env, do_active=False, gain=None, - do_annots=False, do_transits=False, - annotation_file=None, audible_transits_pkl=None, - terraced=False, fill_layers=False, max_tracks=1000, - ): - # 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 - - # 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) - - # 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) - - # 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.show() - - def plot_dem(self, show_scalar_bar=False): - # load DEM - dem = load_DEM(self.project_dir, self.unit, self.site) - 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 - - # 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, y, z, color="white"): - 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): - 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): - 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)"]) - 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): - 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, i=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=300, opacity=1): - 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 plot_annotations(self, annotation_file=None): - # load annotations - if annotation_file is None: - 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}") - annotations = Annotations(annotation_file, only_valid=True) - - if annotations.empty: - print("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") - 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") - - # clip to the area - using a box() plays more nicely with z values than the study area itself - annotations = annotations.to_crs(self.crs) - annotations = annotations.clip(box(*self.study_area.total_bounds)).explode() - - # create and plot segments - audible_actors = [] - inaudible_actors = [] - for _, annot in annotations.iterrows(): - color = "deepskyblue" if annot["audible"] else "red" - polyline = create_polyline_3d(annot["geometry"]) - 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) - - 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=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*{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 setup_z_scale(self): - 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 - ) - +"""CLI entry: python -m nps_active_space.scripts.viz or scripts/viz.py path.""" +from nps_active_space.viz.cli import main 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", 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("--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", 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("--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() - - usy = args.deployment - unit, site, year = usy[:4], usy[4:-4], usy[-4:] - 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 - - Visualizer(unit, site, year, args.environment, - do_active, args.gain, - do_annotations, do_transits, - args.annotation_file, args.transits_pkl, - args.terraced, args.fill_layers, args.max_tracks) \ No newline at end of file + main() diff --git a/nps_active_space/utils/ais/query.py b/nps_active_space/utils/ais/query.py index 8b2fd53d..d44b5d74 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/utils/helpers.py b/nps_active_space/utils/helpers.py index 1599650e..feffaece 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: @@ -161,13 +169,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/nps_active_space/utils/load_tracks.py b/nps_active_space/utils/load_tracks.py new file mode 100644 index 00000000..986526a3 --- /dev/null +++ b/nps_active_space/utils/load_tracks.py @@ -0,0 +1,133 @@ +from pathlib import Path +from typing import NamedTuple + +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.enums import TrackSource +from nps_active_space.utils.helpers import create_overflights_engine, query_adsb, query_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 + + +def _faa_paths() -> tuple[str, str]: + return ( + cfg.read("project", "FAA_Releasable_db"), + cfg.read("project", "FAA_type_corrections"), + ) + + +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"), + 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") + # 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( + engine=engine, + start_date=start_date, + end_date=end_date, + mask=study_area, + ) + 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 | None, +) -> 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") + _apply_site_local_timestamps(tracks, microphone) + return LoadedTracks(tracks, None, None) + + +def load_tracks( + source: TrackSource, + *, + start_date: str, + end_date: str, + study_area: gpd.GeoDataFrame, + microphone: Microphone | None = None, + include_faa_paths: bool = True, +) -> LoadedTracks: + """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, + include_faa_paths=include_faa_paths, + ) + case TrackSource.GPS: + 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 _: + raise ValueError(f"Unknown track source: {source}") diff --git a/nps_active_space/viz/__init__.py b/nps_active_space/viz/__init__.py new file mode 100644 index 00000000..621c6abd --- /dev/null +++ b/nps_active_space/viz/__init__.py @@ -0,0 +1,52 @@ +"""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_track_source_args, + 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_track_source_args", + "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 00000000..e6f50de3 --- /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/assets/waypoint_icon.png b/nps_active_space/viz/assets/waypoint_icon.png new file mode 100644 index 00000000..805c5b7f Binary files /dev/null and b/nps_active_space/viz/assets/waypoint_icon.png differ diff --git a/nps_active_space/viz/cli.py b/nps_active_space/viz/cli.py new file mode 100644 index 00000000..4226dfdd --- /dev/null +++ b/nps_active_space/viz/cli.py @@ -0,0 +1,228 @@ +from __future__ import annotations + +import argparse +import sys +from datetime import datetime +from pathlib import Path + +from nps_active_space.utils.enums import TrackSource +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, + track_source: TrackSource | None = None, + plot_all: bool = False, + annotation_file: str | None = None, + transits_pkl: str | None = None, +) -> 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 + 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() + + 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( + "--track-source", + type=TrackSource, + choices=list(TrackSource), + 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", + 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 tracks to show (annotations, audible transits, or --track-source).", + ) + 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="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). Requires --track-source. 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 + + 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=track_source, + 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, + track_source, + 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 00000000..a89d7778 --- /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 00000000..4ed47775 --- /dev/null +++ b/nps_active_space/viz/geometry.py @@ -0,0 +1,114 @@ +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, + *, + 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/markers.py b/nps_active_space/viz/markers.py new file mode 100644 index 00000000..200c029c --- /dev/null +++ b/nps_active_space/viz/markers.py @@ -0,0 +1,56 @@ +"""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). + + 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), + } + + +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 new file mode 100644 index 00000000..6a99eabd --- /dev/null +++ b/nps_active_space/viz/visualizer.py @@ -0,0 +1,623 @@ +from __future__ import annotations + +import glob +import configparser +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.computation import NMSIM_bbox_utm +from nps_active_space.utils.enums import TrackSource +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.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, + apply_window_icon, + 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" + 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 + + 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, + track_source: TrackSource | None = None, + annotation_file: str | None = None, + audible_transits_pkl: 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, + ) -> 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(title=WINDOW_TITLE) + apply_window_icon(self.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 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() + 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 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: + 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, + line_width=4, + render_lines_as_tubes=True, + ) + 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: + self._status("No audible transits pkl file 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_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. + + 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) + + self._status(f"Querying {source} tracks from {start_date} to {end_date}") + try: + loaded = load_tracks( + source, + start_date=start_date, + end_date=end_date, + study_area=self.study_area, + microphone=microphone, + include_faa_paths=False, + ) + 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: + self._status(f"No {source} tracks loaded.") + return + + track_ids = tracks["track_id"].drop_duplicates() + self._status(f"{len(track_ids)} {source} tracks ({len(tracks)} points)") + if len(track_ids) > self.max_tracks: + 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)] + self._status(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, + include_z=source is not TrackSource.AIS, + ) + if line.is_empty or line.length == 0: + continue + match source: + case TrackSource.AIS: + polyline = self._sea_surface_polyline(line) + case TrackSource.ADSB | TrackSource.GPS: + 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: + self._status(f"No {source} 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=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/ground_truthing/test_segments.py b/tests/ground_truthing/test_segments.py index a6f2448d..3dfa9371 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/utils/ais/test_alignment.py b/tests/utils/ais/test_alignment.py index dbcea43e..11e70b16 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) diff --git a/tests/utils/ais/test_query.py b/tests/utils/ais/test_query.py index 26895bdb..34156c13 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 5b6a71d9..f8c02fc6 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_layered_activespace.py b/tests/utils/test_load_layered_activespace.py new file mode 100644 index 00000000..5e93b34d --- /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/utils/test_load_tracks.py b/tests/utils/test_load_tracks.py new file mode 100644 index 00000000..64bb6ef3 --- /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 new file mode 100644 index 00000000..6a97bdee --- /dev/null +++ b/tests/viz/test_viz.py @@ -0,0 +1,535 @@ +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, 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, + 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_track_source_args, + resolve_viz_plot_flags, + sea_surface_z_profile, + track_points_to_linestring, + utm_orientation_axes_kwargs, + vertex_z_from_coord, +) + + +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 TestResolveVizPlotFlags: + def test_annotation_file_implies_annotations(self): + flags = resolve_viz_plot_flags(annotation_file="/tmp/a.geojson") + 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, None) + + def test_all_enables_standard_layers(self): + flags = resolve_viz_plot_flags(plot_all=True) + assert flags == (True, True, True, None) + + @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) + + +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: + 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)]) + 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.viz.elevation.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.viz.elevation.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 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 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_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" + assert kwargs["ylabel"] == "N" + assert kwargs["zlabel"] == "Z" + assert kwargs["viewport"] == (0, 0, 0.2, 0.2) + + +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 + + 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"): + 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 + + 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 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 nps_active_space.viz.cli import main + + assert scripts_viz.main is main