diff --git a/pyproject.toml b/pyproject.toml index 465db797..86c587c5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -59,6 +59,7 @@ dynamic = ["version"] mission = [ "echopype", "echoshader", + "echoregions", "opencv", ] notebook = [ diff --git a/src/echodataflow/deployment/flow_registry.py b/src/echodataflow/deployment/flow_registry.py index 5d1cf0fb..92c9822d 100644 --- a/src/echodataflow/deployment/flow_registry.py +++ b/src/echodataflow/deployment/flow_registry.py @@ -69,8 +69,23 @@ class FlowRegistration: "update_cache_MVBS": FlowRegistration( entrypoint="echodataflow/flows/flows_viz_cloud.py:flow_update_cache_MVBS", ), + "update_cache_CPS": FlowRegistration( + entrypoint="echodataflow/flows/flows_viz_cloud.py:flow_update_cache_CPS", + description="Update the visualization cache from the latest CPS plot product.", + ), "transect_update": FlowRegistration( entrypoint="echodataflow/flows/flows_transect.py:flow_transect_update", description="Process updates to transect start/end information.", ), -} + "process_CPS": FlowRegistration( + entrypoint="echodataflow/flows/flows_CPS.py:flow_process_CPS", + description="Process CPS acoustic data for new or updated transect segments.", + ), + "simulate_transects": FlowRegistration( + entrypoint=( + "echodataflow/flows/flows_simulation.py:" + "flow_simulate_transects" + ), + description="Simulate realtime transect updates for testing.", + ), +} \ No newline at end of file diff --git a/src/echodataflow/flows/flows_CPS.py b/src/echodataflow/flows/flows_CPS.py new file mode 100644 index 00000000..4c8b9846 --- /dev/null +++ b/src/echodataflow/flows/flows_CPS.py @@ -0,0 +1,1047 @@ +import asyncio +from pathlib import Path + +import dask_image.ndfilters +import echopype as ep +import echoregions as er +import numpy as np +import pandas as pd +import xarray as xr +from prefect import flow, get_client, runtime +from prefect.states import Cancelled +from echodataflow.flows.flows_helper import deployment_already_running +from prefect_dask import DaskTaskRunner + +from echodataflow.utils.processing_ledger import get_completed_sv_files, resolve_database +from echodataflow.tasks.tasks_acoustics import ( + task_compute_NASC_from_masked_Sv, +) + + +def _pick_channel_by_frequency( + ds: xr.Dataset, + freq_hz: float, +) -> str: + idx = int( + np.abs( + ds["frequency_nominal"].values - freq_hz + ).argmin() + ) + return str(ds["channel"].values[idx]) + + +def _dilate_7x7( + da: xr.DataArray, +) -> xr.DataArray: + return xr.DataArray( + dask_image.ndfilters.maximum_filter( + da.data, + size=(1, 7, 7), + ), + dims=da.dims, + coords=da.coords, + ) + + +def _mask_above_seafloor( + ds: xr.Dataset, + bottom_line_path: Path, + channel: str, +) -> xr.DataArray: + lines = er.read_lines_csv( + str(bottom_line_path) + ) + + depth_da = ds["depth"].sel( + channel=channel + ) + + deepest_ping_idx = int( + depth_da.max( + dim="range_sample", + skipna=True, + ) + .argmax(dim="ping_time") + .values + ) + + depth = ( + depth_da + .isel(ping_time=deepest_ping_idx) + .dropna( + dim="range_sample", + how="all", + ) + ) + + valid_length = depth.sizes[ + "range_sample" + ] + + ds_trimmed = ds.isel( + range_sample=slice( + 0, + valid_length, + ) + ) + + sv_target = ds_trimmed["Sv"].sel( + channel=[channel] + ) + + sv_for_regions = xr.DataArray( + sv_target.values, + dims=[ + "channel", + "ping_time", + "depth", + ], + coords={ + "channel": [channel], + "ping_time": sv_target["ping_time"], + "depth": depth.values, + }, + name="Sv", + ) + + bottom_mask, _ = lines.seafloor_mask( + sv_for_regions, + operation="above_below", + method="slinear", + limit_area=None, + limit_direction="both", + ) + + above_mask = ~bottom_mask.astype(bool) + + if "channel" in above_mask.dims: + above_mask = above_mask.squeeze( + "channel", + drop=True, + ) + + above_mask = above_mask.rename( + {"depth": "range_sample"} + ) + + above_mask = above_mask.assign_coords( + range_sample=ds_trimmed[ + "range_sample" + ] + ) + + return above_mask.reindex( + range_sample=ds["range_sample"], + fill_value=False, + ) + + +def _export_nasc_to_echoview_csv( + ds_nasc: xr.Dataset, + output_filepath: Path, + process_id: int = 1928, +) -> None: + df = ( + ds_nasc + .to_dataframe() + .reset_index() + ) + + depth_step = ( + float( + np.nanmedian( + np.diff( + ds_nasc["depth"].values + ) + ) + ) + if ds_nasc.sizes.get( + "depth", + 0, + ) > 1 + else 5.0 + ) + + ping_time = pd.to_datetime( + df["ping_time"] + ) + + out = pd.DataFrame( + { + "Process_ID": process_id, + "Interval": ( + pd.factorize( + df["distance"] + )[0] + + 1 + ), + "Layer": ( + pd.factorize( + df["depth"] + )[0] + + 1 + ), + "Sv_mean": -999.0, + "NASC": df["NASC"].fillna(0.0), + "Height_mean": depth_step, + "Depth_mean": df["depth"], + "Layer_depth_min": ( + df["depth"] + - depth_step / 2 + ), + "Layer_depth_max": ( + df["depth"] + + depth_step / 2 + ), + "Ping_S": 0, + "Ping_E": 0, + "Dist_M": ( + df["distance"] * 1852 + ), + "Date_M": ping_time.dt.strftime( + "%Y%m%d" + ), + "Time_M": ( + ping_time + .dt.strftime( + "%H:%M:%S.%f" + ) + .str[:-2] + ), + "Lat_M": df["latitude"], + "Lon_M": df["longitude"], + "Noise_Sv_1m": -999.0, + "Minimum_Sv_threshold_applied": 1, + "Maximum_Sv_threshold_applied": 0, + "Standard_deviation": 0.0, + "Thickness_mean": depth_step, + "Range_mean": df["depth"], + "Exclude_below_line_range_mean": 999.0, + "Exclude_above_line_range_mean": 0.0, + } + ) + + output_filepath.parent.mkdir( + parents=True, + exist_ok=True, + ) + + out.to_csv( + output_filepath, + index=False, + encoding="utf-8-sig", + ) + + +@flow( + log_prints=True, + task_runner=DaskTaskRunner(), +) +def flow_process_CPS( + path_transect_csv: str, + path_snapshot_csv: str, + path_main: str, + processing_db: str = "processing.db", + target_frequency: float = 70000, + min_depth: float = 10.0, + seafloor_threshold: list = [-40, 2.4, 1.0], + seafloor_offset: float = 0.5, + seafloor_r0: float = 10, + seafloor_r1: float = 1000, + seafloor_wtheta: int = 28, + seafloor_wphi: int = 52, + mask_mode: str = "cps", + fallback_sv_threshold: float = -70, + range_bin: str = "10m", + dist_bin: str = "0.5nmi", + nasc_process_id: int = 1928, +): + + # Prevent overlapping runs of this deployment + already_running = asyncio.run( + deployment_already_running() + ) + + if already_running: + + async def cancel_run(): + async with get_client() as client: + await client.set_flow_run_state( + flow_run_id=runtime.flow_run.id, + state=Cancelled( + message=( + "Another instance of this " + "flow is already running" + ) + ), + ) + + asyncio.run(cancel_run()) + return + + path_main = Path(path_main) + + path_transect = Path( + path_transect_csv + ) + + path_snapshot = Path( + path_snapshot_csv + ) + + path_sv = ( + path_main / "Sv" + ) + + db_path = resolve_database(path_main, processing_db) + + path_cps = ( + path_main / "CPS_Masks_Zarr" + ) + + path_bottom = ( + path_main / "CPS_Seafloor_CSVs" + ) + + path_nasc = ( + path_main / "CPS_NASC_Zarr" + ) + + path_nasc_csv = ( + path_main / "CPS_NASC_CSV" + ) + + for path in [ + path_cps, + path_bottom, + path_nasc, + path_nasc_csv, + ]: + path.mkdir( + parents=True, + exist_ok=True, + ) + + # --------------------------------------------- + # Find completed transects still needing CPS + # --------------------------------------------- + + current = pd.read_csv( + path_transect, + dtype={ + "transectPart": "string", + "transectNumber": "string", + "transectStart": "string", + "transectEnd": "string", + }, + ) + + # Ignore transects that have not finished yet. + completed = current.dropna( + subset=[ + "transectPart", + "transectStart", + "transectEnd", + ] + ).copy() + + pending_rows = [] + + for _, transect in completed.iterrows(): + + name = ( + f"transect_" + f"{transect['transectPart']}" + ) + + cps_output = ( + path_cps + / f"{name}_CPS.zarr" + ) + + nasc_output = ( + path_nasc + / f"{name}_nasc.zarr" + ) + + # A transect is considered complete only when + # both CPS and NASC products exist. + if ( + cps_output.exists() + and nasc_output.exists() + ): + continue + + pending_rows.append( + transect + ) + + changed = pd.DataFrame( + pending_rows, + columns=current.columns, + ) + + if changed.empty: + current.to_csv( + path_snapshot, + index=False, + ) + print( + "No completed transects require CPS processing." + ) + return + + if isinstance(db_path, Path) and not db_path.exists(): + print( + f"Processing ledger not found: " + f"{db_path}" + ) + return + + # --------------------------------------------- + # Process each changed transect + # --------------------------------------------- + + for _, transect in changed.iterrows(): + + start = pd.to_datetime( + transect["transectStart"], + utc=True, + ) + + end = pd.to_datetime( + transect["transectEnd"], + utc=True, + ) + + name = ( + f"transect_" + f"{transect['transectPart']}" + ) + + sv_filenames = get_completed_sv_files( + db_path, + start_time=start, + end_time=end, + ) + + if not sv_filenames: + print( + f"No Sv data for {name}" + ) + continue + + # ----------------------------------------- + # Build continuous Sv transect + # ----------------------------------------- + + sv_paths = [ + path_sv / filename + for filename in sv_filenames + ] + + sv_paths = [ + path + for path in sv_paths + if path.exists() + ] + + if not sv_paths: + continue + + datasets = [ + xr.open_zarr(path) + for path in sv_paths + ] + + ds = xr.concat( + datasets, + dim="ping_time", + data_vars="minimal", + coords="minimal", + compat="override", + ).sortby( + "ping_time" + ) + + _, unique_idx = np.unique( + ds["ping_time"].values, + return_index=True, + ) + + ds = ds.isel( + ping_time=np.sort(unique_idx) + ) + + if ds.sizes.get("ping_time", 0) == 0: + continue + + # ----------------------------------------- + # Require complete Sv coverage + # ----------------------------------------- + + expected_start = start.tz_convert(None) + expected_end = end.tz_convert(None) + + coverage_start = pd.Timestamp( + ds["ping_time"].values[0] + ) + + coverage_end = pd.Timestamp( + ds["ping_time"].values[-1] + ) + + tolerance = pd.Timedelta(seconds=5) + + if ( + coverage_start > expected_start + tolerance + or coverage_end < expected_end - tolerance + ): + print( + f"{name}: incomplete Sv coverage. " + f"Available: {coverage_start} -> {coverage_end}; " + f"required: {expected_start} -> {expected_end}. " + "Leaving transect pending." + ) + continue + + # We know the complete transect is now covered. + ds = ds.sel( + ping_time=slice( + expected_start, + expected_end, + ) + ) + + print( + f"{name}: " + f"{len(sv_paths)} Sv files, " + f"{ds.sizes['ping_time']} pings" + ) + + # ----------------------------------------- + # CPS processing + # ----------------------------------------- + + chunks = { + "channel": 1, + "ping_time": 1000, + "range_sample": -1, + } + + target_channel = ( + _pick_channel_by_frequency( + ds, + target_frequency, + ) + ) + + chunked = ds.chunk( + chunks + ) + + # ----------------------------------------- + # Common geometry + # ----------------------------------------- + + aligned = ( + ep.commongrid + .resample_to_geometry( + chunked, + target_variable="Sv", + target_channel=target_channel, + ) + ) + + if "sound_absorption" in ds: + aligned[ + "sound_absorption" + ] = ds[ + "sound_absorption" + ] + + aligned = ( + ep.consolidate.add_depth( + aligned + ) + ) + + ds[ + ["Sv", "echo_range"] + ] = aligned[ + ["Sv", "echo_range"] + ] + + ds = ( + ep.consolidate.add_depth( + ds + ) + ) + + # ----------------------------------------- + # Background noise + # ----------------------------------------- + + try: + ds = ( + ep.clean + .remove_background_noise( + ds, + ping_num=20, + range_sample_num=5, + SNR_threshold="5.0dB", + ) + ) + + except Exception as exc: + print( + f"{name}: background-noise " + f"removal failed: {exc}" + ) + + ds["Sv_corrected"] = ( + ds["Sv"] + ) + + sv_var = ( + "Sv_corrected" + if "Sv_corrected" in ds + else "Sv" + ) + + # ----------------------------------------- + # Detect seafloor with Blackwell + # ----------------------------------------- + + bottom_path = None + + try: + bottom = ( + ep.mask.detect_seafloor( + ds=ds, + method="blackwell", + params={ + "channel": ( + target_channel + ), + "var_name": "Sv", + "threshold": ( + seafloor_threshold + ), + "offset": ( + seafloor_offset + ), + "r0": ( + seafloor_r0 + ), + "r1": ( + seafloor_r1 + ), + "wtheta": ( + seafloor_wtheta + ), + "wphi": ( + seafloor_wphi + ), + }, + ) + ) + + bottom_df = pd.DataFrame( + { + "time": ( + bottom[ + "ping_time" + ].values + ), + "depth": ( + bottom.values + ), + } + ) + + bottom_df = ( + bottom_df[ + bottom_df[ + "depth" + ] > -0.2 + ] + ) + + bottom_path = ( + path_bottom + / f"{name}_bottom_line.csv" + ) + + bottom_df.to_csv( + bottom_path, + index=False, + ) + + except Exception as exc: + print( + f"{name}: seafloor " + f"detection failed: {exc}" + ) + + # ----------------------------------------- + # Build valid water-column mask + # + # 1. Exclude surface <= min_depth + # 2. Exclude seafloor and everything below + # + # This happens BEFORE CPS classification. + # ----------------------------------------- + + target_depth = ( + ds["depth"].sel( + channel=target_channel + ) + ) + + surface_mask = ( + target_depth > min_depth + ) + + # If seafloor detection/masking fails, + # keep the entire water column apart from + # the surface exclusion. + above_seafloor_mask = ( + xr.ones_like( + surface_mask, + dtype=bool, + ) + ) + + if bottom_path is not None: + try: + above_seafloor_mask = ( + _mask_above_seafloor( + ds, + bottom_path, + target_channel, + ) + ) + + except Exception as exc: + print( + f"{name}: echoregions " + f"seafloor mask failed: {exc}" + ) + + valid_water_column = ( + surface_mask + & above_seafloor_mask + ) + + # Save intermediate masks/products for diagnostics + ds["surface_mask"] = surface_mask + ds["above_seafloor_mask"] = above_seafloor_mask + ds["valid_water_column"] = valid_water_column + + ds["Sv_water_column"] = ( + ds["Sv"].where(valid_water_column) + ) + + # Broadcast the 2-D water-column mask + # (ping_time, range_sample) over channels. + # + # CPS calculations below therefore never + # see the upper 10 m or the seafloor. + sv_for_cps = ( + ds[sv_var].where( + valid_water_column + ) + ) + + # ----------------------------------------- + # CPS classifier + # ----------------------------------------- + + try: + + # Smooth ONLY the valid water column + ds["Sv_smoothed"] = ( + sv_for_cps + .rolling( + ping_time=3, + range_sample=11, + ) + .mean() + ) + + # Variance using the already-masked + # water-column Sv + ds["variance"] = ( + 10 ** ( + sv_for_cps / 10 + ) + - 10 ** ( + ds["Sv_smoothed"] + / 10 + ) + ) ** 2 + + ds["variance_smoothed"] = ( + ds["variance"] + .rolling( + ping_time=3, + range_sample=11, + ) + .mean() + ) + + ds["variance_smoothed"] = ( + 10 + * np.log10( + ds[ + "variance_smoothed" + ] + ** 0.5 + ) + ) + + ds["variance_smoothed"] = ( + _dilate_7x7( + ds[ + "variance_smoothed" + ] + ) + ) + + if ( + mask_mode == "cps" + and ds.sizes[ + "channel" + ] >= 4 + ): + + ch38 = ( + _pick_channel_by_frequency( + ds, + 38000, + ) + ) + + ch70 = ( + _pick_channel_by_frequency( + ds, + 70000, + ) + ) + + ch120 = ( + _pick_channel_by_frequency( + ds, + 120000, + ) + ) + + ch200 = ( + _pick_channel_by_frequency( + ds, + 200000, + ) + ) + + # ----------------------------- + # Variance criteria + # ----------------------------- + + sd200 = ( + ds[ + "variance_smoothed" + ].sel( + channel=ch200 + ) + ) + + sd120 = ( + ds[ + "variance_smoothed" + ].sel( + channel=ch120 + ) + ) + + mask_sd = ( + (sd200 > -65) + & (sd120 > -65) + ) + + # ----------------------------- + # Frequency-response criteria + # + # This is now calculated AFTER + # surface + bottom removal. + # ----------------------------- + + ds["Sv_dilated"] = ( + _dilate_7x7( + ds["Sv_smoothed"] + ) + ) + + diff = ( + ds["Sv_dilated"] + - ds[ + "Sv_dilated" + ].sel( + channel=ch38 + ) + ) + + mask_frequency = ( + ( + diff.sel( + channel=ch200 + ) + > -13.51 + ) + & ( + diff.sel( + channel=ch200 + ) + < 12.53 + ) + & ( + diff.sel( + channel=ch120 + ) + > -13.50 + ) + & ( + diff.sel( + channel=ch120 + ) + < 9.37 + ) + & ( + diff.sel( + channel=ch70 + ) + > -13.85 + ) + & ( + diff.sel( + channel=ch70 + ) + < 9.89 + ) + ) + + final_mask = ( + mask_frequency + & mask_sd + & valid_water_column + ) + + else: + + # Fallback also uses ONLY the + # valid water column. + final_mask = ( + sv_for_cps + > fallback_sv_threshold + ) + + except Exception as exc: + + print( + f"{name}: CPS classifier " + f"failed: {exc}" + ) + + # Same rule for fallback: + # surface and bottom remain excluded. + final_mask = ( + sv_for_cps + > fallback_sv_threshold + ) + + # ----------------------------------------- + # Apply final CPS mask + # + # Values come from original Sv. + # Classification was performed on + # pre-masked / cleaned water-column Sv. + # ----------------------------------------- + + ds["Sv_masked"] = ( + ds["Sv"].where( + final_mask + ) + ) + + # ----------------------------------------- + # Save CPS product + # ----------------------------------------- + + cps_path = ( + path_cps + / f"{name}_CPS.zarr" + ) + + for variable in ds.variables: + ds[ + variable + ].encoding.pop( + "chunks", + None, + ) + + ds.chunk( + chunks + ).to_zarr( + cps_path, + mode="w", + consolidated=True, + ) + + # ----------------------------------------- + # NASC + # + # Sv_masked already contains ONLY: + # + # depth > min_depth + # above seafloor + # CPS-positive samples + # + # Therefore no additional surface or + # bottom masking is necessary here. + # ----------------------------------------- + + ds_nasc = ( + task_compute_NASC_from_masked_Sv( + ds_Sv_masked=ds, + range_bin=range_bin, + dist_bin=dist_bin, + ) + ) + + nasc_path = ( + path_nasc + / f"{name}_nasc.zarr" + ) + + ds_nasc.to_zarr( + nasc_path, + mode="w", + consolidated=True, + ) + + _export_nasc_to_echoview_csv( + ds_nasc, + path_nasc_csv + / f"{name}_nasc.csv", + process_id=nasc_process_id, + ) + + print( + f"{name}: CPS + NASC complete" + ) + + current.to_csv( + path_snapshot, + index=False, + ) \ No newline at end of file diff --git a/src/echodataflow/flows/flows_acoustics.py b/src/echodataflow/flows/flows_acoustics.py index 2659db39..493126df 100644 --- a/src/echodataflow/flows/flows_acoustics.py +++ b/src/echodataflow/flows/flows_acoustics.py @@ -12,6 +12,7 @@ from prefect.futures import as_completed from prefect.states import Cancelled, Failed from prefect import runtime +from prefect.events import emit_event from echodataflow.flows.flows_helper import deployment_already_running from echodataflow.deployment.task_runners import dask_task_runner_from_environment @@ -77,6 +78,9 @@ def flow_raw2Sv( path_main: str = "", processing_db: str = "processing.db", new_file_num_limit: int = 50, + add_depth: bool = True, + add_location: bool = True, + add_splitbeam_angle: bool = False, ): # Check if the deployment is already running already_running = asyncio.run(deployment_already_running()) @@ -144,6 +148,9 @@ async def cancel_run(): sonar_model=sonar_model, datagram_type=datagram_type, nmea_sentence=nmea_sentence, + add_depth=add_depth, + add_location=add_location, + add_splitbeam_angle=add_splitbeam_angle, ) errors = [] @@ -245,6 +252,14 @@ async def set_failed_state(): asyncio.run(set_failed_state()) raise RuntimeError(error_msg) + emit_event( + event="echodataflow.sv.updated", + resource={ + "prefect.resource.id": "sv-monitor", + "prefect.resource.name": "sv-monitor", + }, + ) + @flow(log_prints=True) async def flow_create_MVBS( diff --git a/src/echodataflow/flows/flows_simulation.py b/src/echodataflow/flows/flows_simulation.py index 1d43582b..4c95e8e6 100644 --- a/src/echodataflow/flows/flows_simulation.py +++ b/src/echodataflow/flows/flows_simulation.py @@ -40,6 +40,7 @@ def flow_copy_raw( path_copy: str = "", s3_bucket: str = "noaa-wcsd-pds", exclude_before: str | None = None, + exclude_after: str | None = None, endpoint_url: str = "https://sdsc.osn.xsede.org", ) -> list[S3CopyResult]: """Copy raw files whose timestamps simulate new realtime arrivals.""" @@ -69,15 +70,28 @@ def flow_copy_raw( # Find the last files that would have been generated # between the previous and current flow execution times idx_wanted = df_raw["timestamp"] < flow_time_curr + if exclude_before is not None: exclude_before_datetime = pd.to_datetime(exclude_before, utc=True) - idx_wanted &= df_raw["timestamp"] > exclude_before_datetime + idx_wanted &= df_raw["timestamp"] >= exclude_before_datetime + + if exclude_after is not None: + exclude_after_datetime = pd.to_datetime(exclude_after, utc=True) + idx_wanted &= df_raw["timestamp"] < exclude_after_datetime + if flow_time_prev is not None: idx_wanted &= df_raw["timestamp"] > flow_time_prev df_raw = df_raw[idx_wanted] if df_raw.empty: print("No new files generated since the last flow execution. Skipping file copy.") + + Variable.set( + _var_key(prefix="prev_start_time"), + flow_time_curr.isoformat(), + overwrite=True, + ) + return [] # Setting up task to download @@ -226,7 +240,7 @@ def flow_simulate_transects( start_transect_num: int = 1, max_transects: int = 20, ) -> None: - """Simulate realtime opening and closing of transects.""" + """Simulate realtime arrival of completed transect rows.""" path_transect = Path(path_transect_csv) path_transect.parent.mkdir( @@ -244,15 +258,11 @@ def flow_simulate_transects( default=None, ) - # --------------------------------------------- - # First run: open first transect - # --------------------------------------------- - if state is None: - transect_num_curr = start_transect_num - action = "open" - else: - transect_num_curr, action = state.split(":") - transect_num_curr = int(transect_num_curr) + transect_num_curr = ( + start_transect_num + if state is None + else int(state) + ) if transect_num_curr > max_transects: print("All simulated transects have been generated.") @@ -279,15 +289,25 @@ def flow_simulate_transects( ) if path_transect.exists(): - df = pd.read_csv( - path_transect, - dtype={ - "transectPart": "string", - "transectNumber": "string", - "transectStart": "string", - "transectEnd": "string", - }, - ) + try: + df = pd.read_csv( + path_transect, + dtype={ + "transectPart": "string", + "transectNumber": "string", + "transectStart": "string", + "transectEnd": "string", + }, + ) + except pd.errors.EmptyDataError: + df = pd.DataFrame( + columns=[ + "transectPart", + "transectNumber", + "transectStart", + "transectEnd", + ] + ) else: df = pd.DataFrame( columns=[ @@ -298,62 +318,21 @@ def flow_simulate_transects( ] ) - # --------------------------------------------- - # OPEN transect - # --------------------------------------------- - if action == "open": - row = pd.DataFrame( - [ - { - "transectPart": transect_num, - "transectNumber": transect_num, - "transectStart": transect_start.isoformat(), - "transectEnd": pd.NA, - } - ] - ) - - df = pd.concat( - [df, row], - ignore_index=True, - ) - - df.to_csv( - path_transect, - index=False, - ) - - Variable.set( - _var_key(prefix="transect_state"), - f"{transect_num_curr}:close", - overwrite=True, - ) - - print( - f"Opened simulated transect {transect_num}: " - f"{transect_start}" - ) - - return - - # --------------------------------------------- - # CLOSE transect - # --------------------------------------------- - idx = ( - df["transectPart"] - == transect_num + row = pd.DataFrame( + [ + { + "transectPart": transect_num, + "transectNumber": transect_num, + "transectStart": transect_start.isoformat(), + "transectEnd": transect_end.isoformat(), + } + ] ) - if not idx.any(): - raise ValueError( - f"Cannot close transect {transect_num}: " - "transect not found in CSV." - ) - - df.loc[ - idx, - "transectEnd", - ] = transect_end.isoformat() + df = pd.concat( + [df, row], + ignore_index=True, + ) df.to_csv( path_transect, @@ -361,12 +340,12 @@ def flow_simulate_transects( ) print( - f"Closed simulated transect {transect_num}: " - f"{transect_end}" + f"Wrote simulated transect {transect_num}: " + f"{transect_start} to {transect_end}" ) Variable.set( _var_key(prefix="transect_state"), - f"{transect_num_curr + 1}:open", + str(transect_num_curr + 1), overwrite=True, ) \ No newline at end of file diff --git a/src/echodataflow/flows/flows_transect.py b/src/echodataflow/flows/flows_transect.py index d3c9ad19..06399962 100644 --- a/src/echodataflow/flows/flows_transect.py +++ b/src/echodataflow/flows/flows_transect.py @@ -53,7 +53,15 @@ def flow_transect_update( # Read the current transect information, preserving transect identifiers # as strings so values with leading zeros (e.g., "002") are not converted # to integers by pandas - current = pd.read_csv(path_transect, dtype="string") + try: + current = pd.read_csv( + path_transect, + dtype="string", + ) + except pd.errors.EmptyDataError: + current = pd.DataFrame( + columns=TRANSECT_COLUMNS + ) if not path_snapshot.exists(): print("No previous transect snapshot found. Initializing snapshot.") diff --git a/src/echodataflow/flows/flows_viz_cloud.py b/src/echodataflow/flows/flows_viz_cloud.py index 96642c9a..b9a1eb71 100644 --- a/src/echodataflow/flows/flows_viz_cloud.py +++ b/src/echodataflow/flows/flows_viz_cloud.py @@ -3,6 +3,7 @@ import configparser import pandas as pd +import numpy as np import xarray as xr import s3fs @@ -123,4 +124,250 @@ def flow_update_cache_MVBS( Path(path_cache) / file_MVBS_zarr, # cache is local mode="w", consolidated=True, - ) \ No newline at end of file + ) + +def _prepare_sv_for_echogram( + ds: xr.Dataset, + var_name: str = "Sv_masked", +) -> xr.Dataset: + """Prepare an Sv dataset for Echoshader visualization.""" + + plot_ds = xr.Dataset( + { + "Sv": ds["Sv"], + "Sv_masked": ds[var_name], + } + ) + + if "Sv_water_column" in ds: + plot_ds["Sv_water_column"] = ds["Sv_water_column"] + + if "frequency_nominal" in ds: + plot_ds["frequency_nominal"] = ds["frequency_nominal"] + + if "depth" in ds: + vertical = ds["depth"] + elif "echo_range" in ds: + vertical = ds["echo_range"] + else: + vertical = ds["range_sample"] + + reduce_dims = [ + dim + for dim in ["channel", "ping_time"] + if dim in vertical.dims + ] + + if reduce_dims: + vertical = vertical.median( + dim=reduce_dims, + skipna=True, + ) + + vertical_values = np.asarray(vertical.values) + + # Sv datasets may have trailing range samples with no valid depth. + # Keep the contiguous valid part used for visualization. + finite = np.isfinite(vertical_values) + + if finite.any(): + invalid = np.flatnonzero(~finite) + + if invalid.size: + valid_length = int(invalid[0]) + else: + valid_length = vertical_values.size + + plot_ds = plot_ds.isel( + range_sample=slice(0, valid_length) + ) + + vertical_values = vertical_values[:valid_length] + + if ( + vertical_values.size == 0 + or not np.isfinite(vertical_values).all() + ): + raise ValueError( + "Could not construct a finite vertical coordinate " + "for Sv visualization." + ) + + plot_ds = plot_ds.assign_coords( + echo_range=( + "range_sample", + vertical_values, + ) + ) + + plot_ds = plot_ds.swap_dims( + {"range_sample": "echo_range"} + ) + + vmin = float( + plot_ds["Sv"].min(skipna=True).compute() + ) + vmax = float( + plot_ds["Sv"].max(skipna=True).compute() + ) + + plot_ds["Sv"] = plot_ds["Sv"].assign_attrs( + actual_range=(vmin, vmax) + ) + + masked_vmin = float( + plot_ds["Sv_masked"].min(skipna=True).compute() + ) + masked_vmax = float( + plot_ds["Sv_masked"].max(skipna=True).compute() + ) + + plot_ds["Sv_masked"] = plot_ds["Sv_masked"].assign_attrs( + actual_range=(masked_vmin, masked_vmax) + ) + + for var in plot_ds.variables: + plot_ds[var].encoding.pop("chunks", None) + plot_ds[var].encoding.pop( + "preferred_chunks", + None, + ) + + return plot_ds + +@flow() +def flow_update_cache_CPS( + path_CPS: str, + path_cache: str, + path_transect_csv: str, + file_CPS_zarr: str = "latest_CPS.zarr", +): + """Update visualization cache from the latest CPS transect product.""" + + path_CPS = Path(path_CPS) + path_cache = Path(path_cache) + path_transect_csv = Path(path_transect_csv) + + path_cache.mkdir( + parents=True, + exist_ok=True, + ) + + # ----------------------------------------------------- + # Find latest completed CPS transect + # ----------------------------------------------------- + + cps_files = sorted( + path_CPS.glob("transect_*_CPS.zarr"), + key=lambda path: path.stat().st_mtime, + ) + + if not cps_files: + print( + f"CPS cache not updated: " + f"no CPS transects found in {path_CPS}" + ) + return + + latest_cps = cps_files[-1] + + transect_number = ( + latest_cps.name + .replace("transect_", "") + .replace("_CPS.zarr", "") + ) + + # ----------------------------------------------------- + # Read transect metadata + # ----------------------------------------------------- + + df_transect = pd.read_csv( + path_transect_csv, + dtype={ + "transectPart": str, + "transectNumber": str, + }, + ) + + transect_rows = df_transect[ + df_transect["transectPart"] + == transect_number + ] + + if transect_rows.empty: + print( + f"Transect {transect_number} " + f"not found in {path_transect_csv}" + ) + return + + transect_row = transect_rows.iloc[-1] + + transect_start = pd.to_datetime( + transect_row["transectStart"], + utc=True, + ).tz_convert(None) + + transect_end = pd.to_datetime( + transect_row["transectEnd"], + utc=True, + ).tz_convert(None) + + # ----------------------------------------------------- + # Open already assembled CPS transect + # ----------------------------------------------------- + + ds_CPS = xr.open_zarr( + latest_cps, + consolidated=True, + ) + + ds_CPS = _prepare_sv_for_echogram( + ds_CPS, + var_name="Sv_masked", + ) + + ds_CPS.attrs.update( + { + "transect_number": str( + transect_number + ), + "transect_start": str( + transect_start + ), + "transect_end": str( + transect_end + ), + "source_cps_file": ( + latest_cps.name + ), + } + ) + + # ----------------------------------------------------- + # Save visualization cache + # ----------------------------------------------------- + + output_path = ( + path_cache + / file_CPS_zarr + ) + + print( + f"Saving transect " + f"{transect_number} " + f"to CPS visualization cache: " + f"{output_path}" + ) + + ds_CPS.chunk( + { + "channel": 1, + "ping_time": -1, + "echo_range": 4000, + } + ).to_zarr( + output_path, + mode="w", + consolidated=True, + ) \ No newline at end of file diff --git a/src/echodataflow/operations/operations_acoustics.py b/src/echodataflow/operations/operations_acoustics.py index cfb5a66e..a5351dbd 100644 --- a/src/echodataflow/operations/operations_acoustics.py +++ b/src/echodataflow/operations/operations_acoustics.py @@ -29,6 +29,9 @@ class RawToSvSettings: sonar_model: str = "EK80" datagram_type: str | None = None nmea_sentence: str | None = None + add_depth: bool = True + add_location: bool = True + add_splitbeam_angle: bool = False @dataclass(frozen=True) @@ -58,23 +61,39 @@ def convert_raw_to_Sv( raw_file=raw_path, sonar_model=settings.sonar_model, ) + ds_sv = ep.calibrate.compute_Sv( echodata=echodata, waveform_mode=settings.waveform_mode, encode_mode=settings.encode_mode, ) - ds_sv = ep.consolidate.add_depth( - ds=ds_sv, - depth_offset=settings.depth_offset, - ) - echodata["Platform"] = echodata["Platform"].drop_duplicates("time1") - ds_sv = ep.consolidate.add_location( - ds=ds_sv, - echodata=echodata, - datagram_type=settings.datagram_type, - nmea_sentence=settings.nmea_sentence, - ) + + if settings.add_splitbeam_angle: + ds_sv = ep.consolidate.add_splitbeam_angle( + ds_sv, + echodata, + waveform_mode=settings.waveform_mode, + encode_mode=settings.encode_mode, + to_disk=False, + ) + + if settings.add_depth: + ds_sv = ep.consolidate.add_depth( + ds=ds_sv, + depth_offset=settings.depth_offset, + ) + + if settings.add_location: + echodata["Platform"] = echodata["Platform"].drop_duplicates("time1") + ds_sv = ep.consolidate.add_location( + ds=ds_sv, + echodata=echodata, + datagram_type=settings.datagram_type, + nmea_sentence=settings.nmea_sentence, + ) + output_path = Path(settings.output_directory) / f"{raw_path.stem}_Sv.zarr" + ds_sv.to_zarr( store=output_path, mode="w", @@ -176,3 +195,21 @@ def create_MVBS( first_ping_time=pd.to_datetime(ds_MVBS["ping_time"][0].values), last_ping_time=pd.to_datetime(ds_MVBS["ping_time"][-1].values), ) + + +def compute_NASC_from_masked_Sv( + ds_Sv_masked: xr.Dataset, + range_bin: str = "10m", + dist_bin: str = "0.5nmi", +) -> xr.Dataset: + """Compute NASC from a masked Sv dataset.""" + + ds_for_nasc = ds_Sv_masked.assign( + Sv=ds_Sv_masked["Sv_masked"] + ) + + return ep.commongrid.compute_NASC( + ds_Sv=ds_for_nasc, + range_bin=range_bin, + dist_bin=dist_bin, + ) \ No newline at end of file diff --git a/src/echodataflow/services/viz_all_cps.py b/src/echodataflow/services/viz_all_cps.py new file mode 100644 index 00000000..b786c8db --- /dev/null +++ b/src/echodataflow/services/viz_all_cps.py @@ -0,0 +1,92 @@ +from pathlib import Path + +import echoshader +import matplotlib.pyplot as plt +import panel as pn +import xarray as xr +from holoviews import opts + +pn.config.autoreload = False + +path_MVBS = Path(r"PATH_TO_YOUR_MVBS_ZARR") +path_CPS = Path(r"PATH_TO_YOUR_CPS_ZARR") + + +def multi_freq_app(): + ds = xr.open_zarr(path_MVBS / "latest_MVBS.zarr") + + egram = ds.eshader.echogram( + channel=list(ds.channel.values), + vmin=-70, + vmax=-36, + cmap="viridis", + opts=opts.Image( + width=1000, + height=400, + tools=["pan", "box_zoom", "wheel_zoom", "reset"], + ), + ) + + return pn.pane.HoloViews(egram) + + +def tricolor_app(): + ds = xr.open_zarr(path_MVBS / "latest_MVBS.zarr") + channels = list(ds.channel.values) + + tricolor = ds.eshader.echogram( + channel=[channels[3], channels[1], channels[0]], + vmin=-70, + vmax=-36, + rgb_composite=True, + opts=opts.RGB( + width=1000, + height=400, + tools=["pan", "box_zoom", "wheel_zoom", "reset"], + ), + ) + + return pn.pane.HoloViews(tricolor) + + +def cps_matplotlib_app(): + zarr_path = next(path_CPS.glob("*_mask_plot.zarr")) + ds = xr.open_zarr(zarr_path) + + channel = "WBT 400142-15 ES70-7C_ES" + da = ds["Sv"].sel(channel=channel).compute() + + x = da["ping_time"].values + y = da["echo_range"].values + z = da.values.T + + fig, ax = plt.subplots(figsize=(13, 5)) + im = ax.pcolormesh(x, y, z, shading="auto", vmin=-115, vmax=-36) + ax.invert_yaxis() + ax.set_xlabel("Ping time") + ax.set_ylabel("Depth / range (m)") + ax.set_title(f"CPS masked Sv - {channel}") + fig.colorbar(im, ax=ax, label="Masked Sv (dB re 1 m$^{-1}$)") + fig.tight_layout() + + return pn.Column( + "# CPS masked Sv", + f"File: `{zarr_path.name}` | Channel: `{channel}`", + pn.pane.Matplotlib(fig, tight=True), + ) + + +test_server = pn.serve( + { + "multi_freq_echogram": multi_freq_app, + "tricolor_echogram": tricolor_app, + "cps_matplotlib_echogram": cps_matplotlib_app, + }, + port=1802, + websocket_origin="*", + admin=True, + show=False, + autoreload=False, + keep_alive=40000, + check_unused_sessions_milliseconds=30000, +) \ No newline at end of file diff --git a/src/echodataflow/services/viz_echogram_track_cps.py b/src/echodataflow/services/viz_echogram_track_cps.py new file mode 100644 index 00000000..8b4806d6 --- /dev/null +++ b/src/echodataflow/services/viz_echogram_track_cps.py @@ -0,0 +1,964 @@ +from pathlib import Path + +import datetime +import os +from sqlalchemy import create_engine, inspect +from echodataflow.utils.processing_ledger import resolve_database + +import holoviews as hv +import numpy as np +import pandas as pd +import panel as pn +import xarray as xr +from holoviews.operation.datashader import rasterize + +hv.extension("bokeh") +pn.extension("tabulator") +pn.config.autoreload = False + + +# --------------------------------------------------------------------- +# Paths / settings +# --------------------------------------------------------------------- + +ROOT_ENV = "ECHODATAFLOW_CPS_ROOT" + +if ROOT_ENV not in os.environ: + raise RuntimeError( + f"{ROOT_ENV} is not set. " + "Set it to the root directory containing the CPS workflow outputs." + ) + +ROOT = Path(os.environ[ROOT_ENV]).expanduser().resolve() + +PATH_CACHE = ROOT / "viz_cache_CPS" +PATH_CPS = ROOT / "CPS_Masks_Zarr" +PATH_NASC = ROOT / "CPS_NASC_Zarr" +PATH_BOTTOM = ROOT / "CPS_Seafloor_CSVs" +PROCESSING_DB = os.environ.get( + "ECHODATAFLOW_CPS_PROCESSING_DB", + "processing.db", +) + +PATH_DB = resolve_database( + ROOT, + PROCESSING_DB, +) +PATH_TRANSECTS = ROOT / "plotSurvey_Survey_Data_Visualizer.csv" + +TARGET_FREQUENCY = float( + os.environ.get("ECHODATAFLOW_CPS_TARGET_FREQUENCY", "70000") +) + + +def pick_channel_by_frequency( + ds: xr.Dataset, + freq_hz: float, +) -> str: + """Return the channel whose nominal frequency is closest to freq_hz.""" + + if "channel" not in ds.coords: + raise ValueError( + "Dataset does not contain a channel coordinate." + ) + + channels = ds["channel"].values + + if "frequency_nominal" not in ds: + return str(channels[0]) + + frequencies = np.asarray( + ds["frequency_nominal"].values + ).squeeze() + + if ( + frequencies.ndim != 1 + or frequencies.size != len(channels) + ): + return str(channels[0]) + + finite = np.isfinite(frequencies) + + if not finite.any(): + return str(channels[0]) + + valid_indices = np.flatnonzero(finite) + + idx = valid_indices[ + np.argmin( + np.abs( + frequencies[finite] - freq_hz + ) + ) + ] + + return str(channels[idx]) + + +# --------------------------------------------------------------------- +# Sv plotting +# --------------------------------------------------------------------- + +def plot_sv( + ds: xr.Dataset, + var_name: str, + channel: str, + title: str, + vmin: float = -100, + vmax: float = -30, +): + """Plot Sv using the original irregular ping times.""" + + da = ds[var_name].sel( + channel=channel + ) + + quadmesh = hv.QuadMesh( + ( + ds["ping_time"].values, + ds["echo_range"].values, + da.values.T, + ), + kdims=[ + "ping_time", + "echo_range", + ], + vdims=["Sv"], + ) + + return rasterize( + quadmesh, + width=1000, + height=400, + ).opts( + cmap="viridis", + clim=( + vmin, + vmax, + ), + invert_yaxis=True, + width=1000, + height=400, + tools=[ + "hover", + "pan", + "box_zoom", + "wheel_zoom", + "reset", + ], + title=title, + ) + + +def plot_seafloor( + transect_number: str, +): + """Load detected seafloor line for a transect.""" + + bottom_path = ( + PATH_BOTTOM + / f"transect_{transect_number}_bottom_line.csv" + ) + + if not bottom_path.exists(): + return None + + bottom = pd.read_csv( + bottom_path + ) + + bottom["time"] = pd.to_datetime( + bottom["time"] + ) + + return hv.Curve( + ( + bottom["time"], + bottom["depth"], + ), + kdims=[ + "ping_time", + ], + vdims=[ + "echo_range", + ], + label="Detected seafloor", + ).opts( + color="black", + line_width=2, + ) + + +# --------------------------------------------------------------------- +# Processing database +# --------------------------------------------------------------------- + +def load_database_tables(): + """Load every table currently present in the processing database.""" + + db_value = str(PATH_DB) + + if "://" in db_value: + database_url = db_value + else: + db_path = Path(db_value) + + if not db_path.exists(): + return {} + + database_url = ( + f"sqlite:///{db_path.resolve().as_posix()}" + ) + + engine = create_engine(database_url) + + inspector = inspect(engine) + table_names = inspector.get_table_names() + + tables = {} + + for table_name in table_names: + tables[table_name] = pd.read_sql_table( + table_name, + con=engine, + ) + + return tables + + +# --------------------------------------------------------------------- +# CPS / NASC product summary +# --------------------------------------------------------------------- + +def load_transect_products(): + """Return a table showing which CPS and NASC products exist.""" + + cps = { + p.name + .replace( + "transect_", + "", + ) + .replace( + "_CPS.zarr", + "", + ) + for p in PATH_CPS.glob( + "transect_*_CPS.zarr" + ) + } + + nasc = { + p.name + .replace( + "transect_", + "", + ) + .replace( + "_nasc.zarr", + "", + ) + for p in PATH_NASC.glob( + "transect_*_nasc.zarr" + ) + } + + transects = sorted( + cps | nasc + ) + + rows = [] + + for transect in transects: + rows.append( + { + "transect": transect, + "CPS": ( + "✓" + if transect in cps + else "" + ), + "NASC": ( + "✓" + if transect in nasc + else "" + ), + } + ) + + return pd.DataFrame( + rows + ) + + +# --------------------------------------------------------------------- +# NASC plotting +# --------------------------------------------------------------------- + +def plot_nasc( + ds_nasc: xr.Dataset, + title: str, +): + """Plot NASC as vertical bars along distance.""" + + if "NASC" not in ds_nasc: + return pn.pane.Markdown( + "### NASC variable not found in dataset" + ) + + nasc = ds_nasc["NASC"] + + # Select the channel closest to the configured target frequency. + if "channel" in nasc.dims: + target_channel = ( + pick_channel_by_frequency( + ds_nasc, + TARGET_FREQUENCY, + ) + ) + + nasc = nasc.sel( + channel=target_channel + ) + + if "frequency_nominal" in nasc.dims: + nasc = nasc.isel( + frequency_nominal=0 + ) + + # Remove singleton dimensions + nasc = nasc.squeeze( + drop=True + ) + + print( + "NASC dims:", + nasc.dims, + ) + print( + "NASC shape:", + nasc.shape, + ) + print( + "NASC coords:", + list( + nasc.coords + ), + ) + + # NASC should ultimately be one value per horizontal interval. + # If another dimension remains, integrate/sum over it for plotting. + while nasc.ndim > 1: + dim_to_reduce = ( + nasc.dims[0] + ) + + nasc = nasc.sum( + dim=dim_to_reduce, + skipna=True, + ) + + dim = nasc.dims[0] + + if "distance" in nasc.coords: + x = nasc[ + "distance" + ].values + xlabel = "Distance (nmi)" + else: + x = nasc[ + dim + ].values + xlabel = dim + + curve = hv.Curve( + ( + x, + nasc.values, + ), + kdims=[ + xlabel, + ], + vdims=[ + "NASC", + ], + ) + + return curve.opts( + width=1000, + height=220, + line_width=2, + tools=[ + "hover", + "pan", + "box_zoom", + "wheel_zoom", + "reset", + ], + title=title, + xlabel=xlabel, + ylabel="NASC", + ) + + +# --------------------------------------------------------------------- +# Latest transect plotting +# --------------------------------------------------------------------- + +def load_latest_transect(): + """Load latest cached CPS dataset.""" + + cache_path = ( + PATH_CACHE + / "latest_CPS.zarr" + ) + + if not cache_path.exists(): + raise FileNotFoundError( + f"CPS cache does not exist yet: " + f"{cache_path}" + ) + + ds = xr.open_zarr( + cache_path + ) + + return ( + cache_path, + ds, + ) + + +def build_latest_transect_panel( + vmin: float = -100, + vmax: float = -30, +): + """Build Original Sv + water-column Sv + CPS masked Sv + NASC.""" + + cache_path, ds = ( + load_latest_transect() + ) + + target_channel = ( + pick_channel_by_frequency( + ds, + TARGET_FREQUENCY, + ) + ) + + target_frequency_label = ( + f"{TARGET_FREQUENCY / 1000:g} kHz" + ) + + transect_number = str( + ds.attrs.get( + "transect_number", + "unknown", + ) + ) + + if transect_number != "unknown": + transect_number = ( + transect_number.zfill( + 3 + ) + ) + + cache_time = ( + datetime.datetime + .fromtimestamp( + cache_path + .stat() + .st_mtime + ) + .strftime( + "%Y-%m-%d %H:%M:%S" + ) + ) + + # -------------------------------------------------------------- + # Transect / Sv time coverage + # -------------------------------------------------------------- + + transect_start = "unknown" + transect_end = "unknown" + + path_transects = ( + PATH_TRANSECTS + ) + + if ( + path_transects.exists() + and transect_number != "unknown" + ): + transect_df = ( + pd.read_csv( + path_transects, + dtype={ + "transectPart": "string", + "transectNumber": "string", + "transectStart": "string", + "transectEnd": "string", + }, + ) + ) + + row = transect_df[ + transect_df[ + "transectNumber" + ] + .str.zfill( + 3 + ) + == transect_number + ] + + if not row.empty: + transect_start = ( + row.iloc[0][ + "transectStart" + ] + ) + + transect_end = ( + row.iloc[0][ + "transectEnd" + ] + ) + + ping_start = ( + pd.to_datetime( + ds[ + "ping_time" + ] + .min() + .values + ) + .strftime( + "%Y-%m-%d %H:%M:%S" + ) + ) + + ping_end = ( + pd.to_datetime( + ds[ + "ping_time" + ] + .max() + .values + ) + .strftime( + "%Y-%m-%d %H:%M:%S" + ) + ) + + # -------------------------------------------------------------- + # Metadata + # -------------------------------------------------------------- + + metadata = pn.pane.Markdown( + f""" +## Latest completed CPS transect + +**Transect:** {transect_number} + +**Transect window:** {transect_start} → {transect_end} + +**Sv coverage:** {ping_start} → {ping_end} UTC + +**Cache updated:** {cache_time} + +**Cache:** `{cache_path.name}` +""" + ) + + # -------------------------------------------------------------- + # Seafloor + # -------------------------------------------------------------- + + bottom_curve = ( + plot_seafloor( + transect_number + ) + ) + + # -------------------------------------------------------------- + # Original Sv + # -------------------------------------------------------------- + + original = plot_sv( + ds, + var_name="Sv", + channel=target_channel, + title=( + f"Original Sv - " + f"{target_frequency_label} | " + f"Transect {transect_number}" + ), + vmin=vmin, + vmax=vmax, + ) + + if bottom_curve is not None: + original = ( + original + * bottom_curve + ) + + # -------------------------------------------------------------- + # Water-column masked Sv + # -------------------------------------------------------------- + + if "Sv_water_column" in ds: + water_column = plot_sv( + ds, + var_name="Sv_water_column", + channel=target_channel, + title=( + f"Water-column masked Sv - " + f"{target_frequency_label} | " + f"Transect {transect_number}" + ), + vmin=vmin, + vmax=vmax, + ) + + if bottom_curve is not None: + water_column = ( + water_column + * bottom_curve + ) + + else: + water_column = ( + pn.pane.Alert( + "Sv_water_column is not available " + "in this CPS product.", + alert_type="warning", + ) + ) + + # -------------------------------------------------------------- + # CPS masked Sv + # -------------------------------------------------------------- + + masked = plot_sv( + ds, + var_name="Sv_masked", + channel=target_channel, + title=( + f"CPS masked Sv - " + f"{target_frequency_label} | " + f"Transect {transect_number}" + ), + vmin=vmin, + vmax=vmax, + ) + + if bottom_curve is not None: + masked = ( + masked + * bottom_curve + ) + + # -------------------------------------------------------------- + # NASC + # -------------------------------------------------------------- + + nasc_path = ( + PATH_NASC + / f"transect_{transect_number}_nasc.zarr" + ) + + if nasc_path.exists(): + try: + ds_nasc = xr.open_zarr( + nasc_path + ) + + nasc_plot = plot_nasc( + ds_nasc, + title=( + f"NASC | " + f"Transect {transect_number}" + ), + ) + + except Exception as e: + nasc_plot = ( + pn.pane.Alert( + f"Could not plot NASC: {e}", + alert_type="warning", + ) + ) + + else: + nasc_plot = ( + pn.pane.Alert( + f"NASC is not available yet for " + f"transect {transect_number}.", + alert_type="info", + ) + ) + + return pn.Column( + metadata, + original, + water_column, + masked, + nasc_plot, + sizing_mode="stretch_width", + ) + + +# --------------------------------------------------------------------- +# Status dashboard +# --------------------------------------------------------------------- + +def build_status_panel(): + """Build live processing status tables.""" + + # -------------------------------------------------------------- + # Product status + # -------------------------------------------------------------- + + product_df = ( + load_transect_products() + ) + + product_table = ( + pn.widgets.Tabulator( + product_df, + pagination=None, + show_index=False, + disabled=True, + sizing_mode="stretch_width", + height=250, + ) + ) + + product_section = ( + pn.Column( + "## Transect products", + product_table, + ) + ) + + # -------------------------------------------------------------- + # processing.db + # -------------------------------------------------------------- + + db_tables = ( + load_database_tables() + ) + + db_panels = [] + + if not db_tables: + db_panels.append( + pn.pane.Alert( + "Processing database is unavailable " + "or contains no tables.", + alert_type="warning", + ) + ) + + else: + for ( + table_name, + df, + ) in db_tables.items(): + table = ( + pn.widgets.Tabulator( + df, + pagination="local", + page_size=10, + show_index=False, + disabled=True, + sizing_mode="stretch_width", + height=300, + ) + ) + + db_panels.append( + pn.Column( + f"### {table_name}", + table, + ) + ) + + database_section = ( + pn.Column( + "## Processing database", + *db_panels, + ) + ) + + update_time = ( + datetime.datetime.now() + .strftime( + "%Y-%m-%d %H:%M:%S" + ) + ) + + header = pn.pane.Markdown( + f""" +# CPS processing monitor + +**Dashboard refreshed:** {update_time} +""" + ) + + return pn.Column( + header, + product_section, + pn.layout.Divider(), + database_section, + sizing_mode="stretch_width", + ) + + +# --------------------------------------------------------------------- +# Main application +# --------------------------------------------------------------------- + +def cps_app(): + """Live CPS monitoring dashboard.""" + + sv_clim = ( + pn.widgets.RangeSlider( + name="Sv color range (dB)", + start=-120, + end=-20, + value=( + -100, + -30, + ), + step=1, + width=450, + ) + ) + + latest_container = ( + pn.Column( + sizing_mode="stretch_width", + ) + ) + + status_container = ( + pn.Column( + sizing_mode="stretch_width", + ) + ) + + def refresh_latest(): + try: + vmin, vmax = ( + sv_clim.value + ) + + latest_container[:] = [ + build_latest_transect_panel( + vmin=vmin, + vmax=vmax, + ) + ] + + print( + "Latest transect panel refreshed at " + f"{datetime.datetime.now():%H:%M:%S}" + ) + + except Exception as e: + latest_container[:] = [ + pn.pane.Alert( + f"Could not load latest CPS " + f"transect: {e}", + alert_type="warning", + ) + ] + + def refresh_status(): + try: + status_container[:] = [ + build_status_panel() + ] + + print( + "Status panel refreshed at " + f"{datetime.datetime.now():%H:%M:%S}" + ) + + except Exception as e: + status_container[:] = [ + pn.pane.Alert( + f"Could not refresh processing " + f"status: {e}", + alert_type="danger", + ) + ] + + # Initial load + refresh_latest() + refresh_status() + + # Rebuild Sv plots when color range changes + sv_clim.param.watch( + lambda event: refresh_latest(), + "value", + ) + + # Refresh echograms / NASC every 30 seconds + pn.state.add_periodic_callback( + refresh_latest, + period=30 * 1000, + ) + + # Refresh processing status every 10 seconds + pn.state.add_periodic_callback( + refresh_status, + period=10 * 1000, + ) + + tabs = pn.Tabs( + ( + "Latest transect", + latest_container, + ), + ( + "Processing status", + status_container, + ), + dynamic=False, + sizing_mode="stretch_width", + ) + + template = ( + pn.template.FastListTemplate( + title="CPS Near-Real-Time Monitor", + main=[ + sv_clim, + tabs, + ], + ) + ) + + return template + + +# --------------------------------------------------------------------- +# Server +# --------------------------------------------------------------------- + +test_server = pn.serve( + { + "cps_echogram": cps_app, + }, + port=1803, + websocket_origin="*", + admin=True, + show=False, + autoreload=False, + keep_alive=40000, + check_unused_sessions_milliseconds=30000, +) \ No newline at end of file diff --git a/src/echodataflow/tasks/tasks_acoustics.py b/src/echodataflow/tasks/tasks_acoustics.py index 153be76b..9e149eb3 100644 --- a/src/echodataflow/tasks/tasks_acoustics.py +++ b/src/echodataflow/tasks/tasks_acoustics.py @@ -1,5 +1,6 @@ """Reusable Prefect tasks for acoustic data processing.""" +import xarray as xr from prefect import get_run_logger, task from echodataflow.operations.operations_acoustics import ( @@ -9,6 +10,7 @@ RawToSvResult, RawToSvSettings, RawToSvWorkItem, + compute_NASC_from_masked_Sv, create_MVBS, convert_raw_to_Sv, ) @@ -32,3 +34,18 @@ def task_create_MVBS( logger = get_run_logger() logger.info(f"Saving MVBS to {item.mvbs_filename}") return create_MVBS(item, settings) + + +@task(log_prints=True) +def task_compute_NASC_from_masked_Sv( + ds_Sv_masked: xr.Dataset, + range_bin: str = "10m", + dist_bin: str = "0.5nmi", +) -> xr.Dataset: + """Compute NASC from a masked Sv dataset as a Prefect task.""" + + return compute_NASC_from_masked_Sv( + ds_Sv_masked=ds_Sv_masked, + range_bin=range_bin, + dist_bin=dist_bin, + ) \ No newline at end of file diff --git a/tests/deployment/test_flow_registry.py b/tests/deployment/test_flow_registry.py index 7bf3ea02..d4cc121a 100644 --- a/tests/deployment/test_flow_registry.py +++ b/tests/deployment/test_flow_registry.py @@ -61,6 +61,15 @@ def test_postprocessing_flows_are_registered_with_realtime_modules(): ) +def test_flow_registry_contains_process_cps_entrypoint(): + registry = importlib.import_module("echodataflow.deployment.flow_registry") + + registration = registry.FLOW_REGISTRY["process_CPS"] + + assert registration.entrypoint == ( + "echodataflow/flows/flows_CPS.py:flow_process_CPS" + ) + def test_resolve_registered_flows_defaults_to_recipe_key( monkeypatch, install_prefect_stubs, diff --git a/tests/test_flow_CPS.py b/tests/test_flow_CPS.py new file mode 100644 index 00000000..b3833baf --- /dev/null +++ b/tests/test_flow_CPS.py @@ -0,0 +1,77 @@ +import pandas as pd + +from echodataflow.flows import flows_CPS + + +async def _false_async(): + return False + + +def test_process_cps_retries_completed_transect_missing_outputs( + monkeypatch, + tmp_path, + capsys, +): + path_main = tmp_path / "output" + path_main.mkdir() + + transect_csv = tmp_path / "transects.csv" + snapshot_csv = tmp_path / "snapshot.csv" + + transect = pd.DataFrame( + { + "transectPart": ["001"], + "transectNumber": ["001"], + "transectStart": ["2024-07-07T00:30:00Z"], + "transectEnd": ["2024-07-07T00:35:00Z"], + } + ) + + # The transect is already present in the snapshot. + # Previously this meant CPS would never reconsider it. + transect.to_csv(transect_csv, index=False) + transect.to_csv(snapshot_csv, index=False) + + # flow_process_CPS requires the processing ledger to exist. + (path_main / "processing.db").touch() + + monkeypatch.setattr( + flows_CPS, + "deployment_already_running", + lambda: _false_async(), + ) + + calls = [] + + def fake_get_completed_sv_files( + db_path, + start_time=None, + end_time=None, + ): + calls.append( + ( + db_path, + start_time, + end_time, + ) + ) + return [] + + monkeypatch.setattr( + flows_CPS, + "get_completed_sv_files", + fake_get_completed_sv_files, + ) + + flows_CPS.flow_process_CPS.fn( + path_transect_csv=str(transect_csv), + path_snapshot_csv=str(snapshot_csv), + path_main=str(path_main), + ) + + output = capsys.readouterr().out + + # Even though the snapshot already contains transect 001, + # it must still be checked because CPS/NASC outputs are missing. + assert len(calls) == 1 + assert "No Sv data for transect_001" in output \ No newline at end of file diff --git a/tests/test_flow_copy_raw.py b/tests/test_flow_copy_raw.py index e55d0e7e..f917ed17 100644 --- a/tests/test_flow_copy_raw.py +++ b/tests/test_flow_copy_raw.py @@ -118,4 +118,45 @@ def test_flow_copy_raw_simulates_new_file_arrivals(monkeypatch, tmp_path): assert any( key.startswith("prev_start_time_") for key in FakeVariable.stored + ) + + +def test_flow_copy_raw_updates_watermark_when_no_files_selected( + monkeypatch, + tmp_path, +): + FakeVariable.stored = {} + + raw_list = tmp_path / "raw_files.csv" + pd.DataFrame( + { + "timestamp": [ + "2024-07-07T00:40:00Z", + ], + "s3_path": [ + "survey/future.raw", + ], + } + ).to_csv(raw_list, index=False) + + monkeypatch.setattr(flows_simulation, "Variable", FakeVariable) + monkeypatch.setattr(flows_simulation.datetime, "datetime", FakeDateTime) + + results = flows_simulation.flow_copy_raw.fn( + path_raw_list=str(raw_list), + path_copy=str(tmp_path / "raw"), + s3_bucket="raw-bucket", + ) + + assert results == [] + + watermark_keys = [ + key + for key in FakeVariable.stored + if key.startswith("prev_start_time_") + ] + + assert len(watermark_keys) == 1 + assert FakeVariable.stored[watermark_keys[0]] == ( + "2024-07-07T00:30:00+00:00" ) \ No newline at end of file diff --git a/tests/test_flow_simulate_transects.py b/tests/test_flow_simulate_transects.py index 87718dd4..d10ff576 100644 --- a/tests/test_flow_simulate_transects.py +++ b/tests/test_flow_simulate_transects.py @@ -16,7 +16,10 @@ def set(cls, key, value, overwrite=False): cls.stored[key] = value -def test_flow_simulate_transects_opens_closes_and_advances(monkeypatch, tmp_path): +def test_flow_simulate_transects_writes_complete_rows_and_advances( + monkeypatch, + tmp_path, +): FakeVariable.stored = {} monkeypatch.setattr(flows_simulation, "Variable", FakeVariable) @@ -30,26 +33,23 @@ def test_flow_simulate_transects_opens_closes_and_advances(monkeypatch, tmp_path "max_transects": 2, } - # First run: open transect 001. + # First run: write complete transect 001. flows_simulation.flow_simulate_transects.fn(**kwargs) df = pd.read_csv(transect_csv, dtype="string") - assert df["transectPart"].tolist() == ["001"] - assert pd.isna(df.loc[0, "transectEnd"]) - - # Second run: close transect 001. - flows_simulation.flow_simulate_transects.fn(**kwargs) - df = pd.read_csv(transect_csv, dtype="string") + assert df["transectPart"].tolist() == ["001"] + assert df.loc[0, "transectStart"] == "2024-07-07T00:00:00+00:00" assert df.loc[0, "transectEnd"] == "2024-07-07T00:10:00+00:00" - # Third run: open transect 002. + # Second run: append complete transect 002. flows_simulation.flow_simulate_transects.fn(**kwargs) df = pd.read_csv(transect_csv, dtype="string") + assert df["transectPart"].tolist() == ["001", "002"] assert df.loc[1, "transectStart"] == "2024-07-07T00:10:00+00:00" - assert pd.isna(df.loc[1, "transectEnd"]) + assert df.loc[1, "transectEnd"] == "2024-07-07T00:20:00+00:00" def test_flow_simulate_transects_stops_after_maximum(monkeypatch, tmp_path, capsys): diff --git a/tests/test_flow_transect.py b/tests/test_flow_transect.py index d7a23921..e3913772 100644 --- a/tests/test_flow_transect.py +++ b/tests/test_flow_transect.py @@ -206,4 +206,72 @@ def test_flow_transect_update_ignores_open_transect(tmp_path, capsys): output = capsys.readouterr().out - assert "No new or updated transect segments." in output \ No newline at end of file + assert "No new or updated transect segments." in output + +def test_get_changed_transects_detects_open_to_closed_update(): + previous = pd.DataFrame( + { + "transectPart": ["001"], + "transectNumber": ["001"], + "transectStart": ["2024-07-07T00:00:00Z"], + "transectEnd": [pd.NA], + } + ) + + current = pd.DataFrame( + { + "transectPart": ["001"], + "transectNumber": ["001"], + "transectStart": ["2024-07-07T00:00:00Z"], + "transectEnd": ["2024-07-07T00:10:00Z"], + } + ) + + changed = get_changed_transects(current, previous) + + assert len(changed) == 1 + assert changed.iloc[0]["transectPart"] == "001" + assert changed.iloc[0]["transectEnd"] == "2024-07-07T00:10:00Z" + +def test_flow_transect_update_handles_header_only_csv( + tmp_path, + capsys, +): + transect_csv = tmp_path / "transects.csv" + snapshot_csv = tmp_path / "snapshot.csv" + path_main = tmp_path / "output" + path_main.mkdir() + + transect_csv.write_text( + "transectPart,transectNumber,transectStart,transectEnd\n" + ) + + flow_transect_update.fn( + path_transect_csv=str(transect_csv), + path_snapshot_csv=str(snapshot_csv), + path_main=str(path_main), + ) + + output = capsys.readouterr().out + + assert snapshot_csv.exists() + assert "No previous transect snapshot found. Initializing snapshot." in output + + +def test_flow_transect_update_handles_zero_byte_csv( + tmp_path, +): + transect_csv = tmp_path / "transects.csv" + snapshot_csv = tmp_path / "snapshot.csv" + path_main = tmp_path / "output" + path_main.mkdir() + + transect_csv.touch() + + flow_transect_update.fn( + path_transect_csv=str(transect_csv), + path_snapshot_csv=str(snapshot_csv), + path_main=str(path_main), + ) + + assert snapshot_csv.exists() \ No newline at end of file