diff --git a/.ci_helpers/docker/docker-compose.yaml b/.ci_helpers/docker/docker-compose.yaml index 0c6130372..8671e7fe0 100644 --- a/.ci_helpers/docker/docker-compose.yaml +++ b/.ci_helpers/docker/docker-compose.yaml @@ -10,12 +10,7 @@ services: volumes: - echopype_miniodata:/data httpserver: - # NOTE ==================================== - # For httpserver test data, - # it needs to be copied to the docker container - # docker cp -L docker_httpserver_1:/usr/local/apache2/htdocs/data ./echopype/test_data - # ===================================== - image: cormorack/http:latest + image: httpd:2.4 ports: - 8080:80 networks: diff --git a/.ci_helpers/docker/http.dockerfile b/.ci_helpers/docker/http.dockerfile deleted file mode 100644 index 96d02eb0c..000000000 --- a/.ci_helpers/docker/http.dockerfile +++ /dev/null @@ -1,4 +0,0 @@ -FROM httpd:2.4 -ARG TARGETPLATFORM - -COPY echopype/test_data /usr/local/apache2/htdocs/data diff --git a/.ci_helpers/docker/setup-services.py b/.ci_helpers/docker/setup-services.py index 559d9bd11..b460e6766 100755 --- a/.ci_helpers/docker/setup-services.py +++ b/.ci_helpers/docker/setup-services.py @@ -6,13 +6,14 @@ import argparse import logging -import shutil +import os import subprocess import sys from pathlib import Path from typing import Dict, List import fsspec +import pooch logger = logging.getLogger("setup-services") streamHandler = logging.StreamHandler(sys.stdout) @@ -24,7 +25,18 @@ HERE = Path(".").absolute() BASE = Path(__file__).parent.absolute() COMPOSE_FILE = BASE / "docker-compose.yaml" -TEST_DATA_PATH = HERE / "echopype" / "test_data" + + +def get_pooch_data_path() -> Path: + """Return path to the Pooch test data cache.""" + ver = os.getenv("ECHOPYPE_DATA_VERSION", "v0.11.1a2") + cache_dir = Path(pooch.os_cache("echopype")) / ver + if not cache_dir.exists(): + raise FileNotFoundError( + f"Pooch cache directory not found: {cache_dir}\n" + "Make sure test data was fetched via conftest.py" + ) + return cache_dir def parse_args(): @@ -77,36 +89,70 @@ def run_commands(commands: List[Dict]) -> None: def load_s3(*args, **kwargs) -> None: + """Populate MinIO with test data from the Pooch cache (skip .zip files).""" + pooch_path = get_pooch_data_path() common_storage_options = dict( client_kwargs=dict(endpoint_url="http://localhost:9000/"), key="minioadmin", secret="minioadmin", ) - bucket_name = "ooi-raw-data" - fs = fsspec.filesystem( - "s3", - **common_storage_options, - ) + bucket_name = "echo-test-data" + fs = fsspec.filesystem("s3", **common_storage_options) test_data = "data" + if not fs.exists(test_data): fs.mkdir(test_data) - if not fs.exists(bucket_name): fs.mkdir(bucket_name) - # Load test data into bucket - for d in TEST_DATA_PATH.iterdir(): - source_path = f"echopype/test_data/{d.name}" - fs.put(source_path, f"{test_data}/{d.name}", recursive=True) + for d in pooch_path.iterdir(): + if d.suffix == ".zip": # skip zip archives to cut redundant I/O + continue + source_path = str(d) + target_path = f"{test_data}/{d.name}" + logger.info(f"Copying {source_path} → {target_path}") + fs.put(source_path, target_path, recursive=True) + + +def load_http_server(http_server_id) -> None: + """Copy test data from Pooch cache to HTTP server container.""" + pooch_path = get_pooch_data_path() + + # Create the data directory in the container + mkdir_result = subprocess.run( + ["docker", "exec", http_server_id, "sh", "-c", "mkdir -p /usr/local/apache2/htdocs/data"], + capture_output=True, + text=True, + ) + if mkdir_result.returncode == 0: + logger.info("Created /usr/local/apache2/htdocs/data directory") + else: + logger.warning(f"mkdir warning (may be harmless): {mkdir_result.stderr}") + + # Copy all dataset directories directly to /usr/local/apache2/htdocs/data/ + for d in pooch_path.iterdir(): + if d.suffix == ".zip": # skip zip archives + continue + source_path = str(d) + dataset_name = d.name + + # Copy directly to the data directory + target_path = f"/usr/local/apache2/htdocs/data/{dataset_name}" + cmd = ["docker", "cp", source_path, f"{http_server_id}:{target_path}"] + logger.info(f"Copying {source_path} → {target_path}") + result = subprocess.run(cmd, capture_output=True, text=True) + + if result.returncode != 0: + logger.error(f"Failed to copy {dataset_name}: {result.stderr}") if __name__ == "__main__": args = parse_args() commands = [] + if all([args.deploy, args.tear_down]): print("Cannot have both --deploy and --tear-down. Exiting.") sys.exit(1) - if not any([args.deploy, args.tear_down]): print("Please provide either --deploy or --tear-down flags. For more help use --help flag.") sys.exit(0) @@ -136,33 +182,29 @@ def load_s3(*args, **kwargs) -> None: } ) - if TEST_DATA_PATH.exists(): - commands.append( - { - "msg": f"Deleting old test folder at {TEST_DATA_PATH} ...", - "cmd": shutil.rmtree, - "args": TEST_DATA_PATH, - } - ) + pooch_path = get_pooch_data_path() + commands.append({"msg": f"Using Pooch test data at {pooch_path}", "cmd": None}) + commands.append( { - "msg": "Copying new test folder from http service ...", - "cmd": [ - "docker", - "cp", - "-L", - f"{args.http_server}:/usr/local/apache2/htdocs/data", - TEST_DATA_PATH, - ], + "msg": "Setting up MinIO S3 bucket with Pooch test data ...", + "cmd": load_s3, } ) - commands.append({"msg": "Setting up minio s3 bucket ...", "cmd": load_s3}) + commands.append( + { + "msg": f"Setting up HTTP server {args.http_server} with Pooch test data ...", + "cmd": load_http_server, + "args": args.http_server, + } + ) if args.tear_down: command = ["docker-compose", "-f", COMPOSE_FILE, "down", "--remove-orphans", "--volumes"] if args.images: - command = command + ["--rmi", "all"] + command += ["--rmi", "all"] commands.append({"msg": "Stopping test services deployment ...", "cmd": command}) + commands.append({"msg": "Done.", "cmd": ["docker", "ps", "--last", "2"]}) run_commands(commands) diff --git a/.ci_helpers/py3.10.yaml b/.ci_helpers/py3.13.yaml similarity index 86% rename from .ci_helpers/py3.10.yaml rename to .ci_helpers/py3.13.yaml index d33ae3ec4..ef2deba36 100644 --- a/.ci_helpers/py3.10.yaml +++ b/.ci_helpers/py3.13.yaml @@ -2,7 +2,7 @@ name: echopype channels: - conda-forge dependencies: - - python=3.10 + - python=3.13 - pip - pip: - -r ../requirements.txt diff --git a/.ci_helpers/setup-services-windows.py b/.ci_helpers/setup-services-windows.py new file mode 100644 index 000000000..dbcab332c --- /dev/null +++ b/.ci_helpers/setup-services-windows.py @@ -0,0 +1,195 @@ +# .ci_helpers/setup-services-windows.py +""" +Spin up local services for Windows CI without Docker: +- Start a MinIO server on localhost:9000 +- Seed S3 from the Pooch cache (unzipped test bundles) +- Start a simple HTTP server on :8080 that exposes ./data (mirrors Linux job) +- Write PID files for clean teardown + +Usage: + python .ci_helpers/setup-services-windows.py start [--no-http] + python .ci_helpers/setup-services-windows.py stop +""" + +import argparse +import os +import pathlib +import shutil +import subprocess +import sys +import time +from urllib.request import urlretrieve + +import fsspec +import pooch + +MINIO_URL = "https://dl.min.io/server/minio/release/windows-amd64/minio.exe" +MINIO_BIN = pathlib.Path(".ci_helpers") / "minio.exe" +STATE_DIR = pathlib.Path(".ci_helpers") / ".state" +STATE_DIR.mkdir(parents=True, exist_ok=True) +MINIO_PID = STATE_DIR / "minio.pid" +HTTP_PID = STATE_DIR / "http.pid" + +# Use localhost everywhere to match tests (which hit http://localhost:9000/...) +MINIO_ENDPOINT = "http://localhost:9000/" +MINIO_USER = "minioadmin" +MINIO_PASS = "minioadmin" + + +def get_pooch_cache() -> pathlib.Path: + """Return the Pooch cache dir for the configured dataset version.""" + ver = os.getenv("ECHOPYPE_DATA_VERSION", "v0.11.1a2") + root = pathlib.Path(pooch.os_cache("echopype")) + path = root / ver + path.mkdir(parents=True, exist_ok=True) + return path + + +def ensure_minio_downloaded() -> None: + if MINIO_BIN.exists(): + return + MINIO_BIN.parent.mkdir(parents=True, exist_ok=True) + print(f"Downloading MinIO -> {MINIO_BIN}", flush=True) + urlretrieve(MINIO_URL, MINIO_BIN) + MINIO_BIN.chmod(0o755) + + +def start_minio() -> None: + """Start MinIO on localhost:9000 and wait until ready.""" + ensure_minio_downloaded() + data_dir = pathlib.Path(os.getenv("USERPROFILE", str(pathlib.Path.home()))) / "minio" / "data" + data_dir.mkdir(parents=True, exist_ok=True) + + env = os.environ.copy() + env["MINIO_ROOT_USER"] = MINIO_USER + env["MINIO_ROOT_PASSWORD"] = MINIO_PASS + + print(f"Starting MinIO on {MINIO_ENDPOINT} (data: {data_dir})", flush=True) + proc = subprocess.Popen( + [ + str(MINIO_BIN), + "server", + str(data_dir), + "--address=localhost:9000", + "--console-address=localhost:9001", + ], + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL, + env=env, + ) + MINIO_PID.write_text(str(proc.pid)) + + # Wait for readiness + import urllib.request + + for _ in range(60): + try: + with urllib.request.urlopen(MINIO_ENDPOINT + "minio/health/ready", timeout=1): + break + except Exception: + time.sleep(1) + else: + raise RuntimeError("MinIO did not become ready on :9000") + + +def seed_s3_from_pooch() -> None: + """Upload unzipped Pooch bundles into MinIO under the 'data/' prefix, + and also mirror a local ./data/ folder for the HTTP server.""" + cache = get_pooch_cache() + + # Seed S3 (MinIO) + fs = fsspec.filesystem( + "s3", + client_kwargs=dict(endpoint_url=MINIO_ENDPOINT), + key=MINIO_USER, + secret=MINIO_PASS, + ) + + for base in ("data", "ooi-raw-data"): + try: + fs.mkdir(base) + except Exception: + pass + + for d in cache.iterdir(): + if d.suffix == ".zip": + continue + tgt = f"data/{d.name}" + print(f"Uploading {d} -> {tgt}", flush=True) + fs.put(str(d), tgt, recursive=True) + + # Build local ./data for HTTP tests that hit http://localhost:8080/data/... + local_data = pathlib.Path("data") + local_data.mkdir(exist_ok=True) + for d in cache.iterdir(): + if d.suffix == ".zip": + continue + dst = local_data / d.name + if not dst.exists(): + if d.is_dir(): + shutil.copytree(d, dst) + else: + shutil.copy2(d, dst) + + +def start_http_server() -> None: + """Serve the repository root so /data/... exists on :8080.""" + root = pathlib.Path(".").absolute() + print(f"Starting local HTTP server on :8080 (root={root})", flush=True) + proc = subprocess.Popen( + [sys.executable, "-m", "http.server", "8080", "--directory", str(root)], + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL, + ) + HTTP_PID.write_text(str(proc.pid)) + + +def stop_pid_file(path: pathlib.Path) -> None: + try: + pid = int(path.read_text().strip()) + except Exception: + return + try: + if sys.platform.startswith("win"): + subprocess.run(["taskkill", "/PID", str(pid), "/F"], check=False) + else: + os.kill(pid, 9) + except Exception: + pass + try: + path.unlink(missing_ok=True) + except Exception: + pass + + +def cmd_start(no_http: bool) -> None: + # Ensure cache exists (prefetch step should have populated it) + _ = get_pooch_cache() + start_minio() + seed_s3_from_pooch() + if not no_http: + start_http_server() + print("Services up.", flush=True) + + +def cmd_stop() -> None: + stop_pid_file(HTTP_PID) + stop_pid_file(MINIO_PID) + print("Services stopped.", flush=True) + + +def main() -> None: + ap = argparse.ArgumentParser() + sub = ap.add_subparsers(dest="cmd", required=True) + s = sub.add_parser("start") + s.add_argument("--no-http", action="store_true", help="Do not start the local HTTP server") + sub.add_parser("stop") + args = ap.parse_args() + if args.cmd == "start": + cmd_start(no_http=args.no_http) + else: + cmd_stop() + + +if __name__ == "__main__": + main() diff --git a/.github/actions/gdrive-rclone/Dockerfile b/.github/actions/gdrive-rclone/Dockerfile deleted file mode 100644 index 948d3ef20..000000000 --- a/.github/actions/gdrive-rclone/Dockerfile +++ /dev/null @@ -1,4 +0,0 @@ -FROM rclone/rclone:latest -RUN apk add bash jq -COPY entrypoint.sh /opt/entrypoint.sh -ENTRYPOINT ["/opt/entrypoint.sh"] diff --git a/.github/actions/gdrive-rclone/action.yaml b/.github/actions/gdrive-rclone/action.yaml deleted file mode 100644 index 8c4e1fde8..000000000 --- a/.github/actions/gdrive-rclone/action.yaml +++ /dev/null @@ -1,5 +0,0 @@ -name: 'GDrive Rclone' -description: 'Performs RClone from a google drive folder' -runs: - using: 'docker' - image: 'Dockerfile' diff --git a/.github/actions/gdrive-rclone/entrypoint.sh b/.github/actions/gdrive-rclone/entrypoint.sh deleted file mode 100755 index 4fb607d5f..000000000 --- a/.github/actions/gdrive-rclone/entrypoint.sh +++ /dev/null @@ -1,35 +0,0 @@ -#!/bin/bash - -set -e - -if [ -d "/opt/rclone_temp" ] -then - echo "/opt/rclone_temp found." -else - echo "creating /opt/rclone_temp" - mkdir /opt/rclone_temp -fi -export RCLONE_DRIVE_SERVICE_ACCOUNT_FILE="/opt/rclone_temp/google-echopype.json" -export RCLONE_DRIVE_ROOT_FOLDER_ID=${ROOT_FOLDER_ID} -export RCLONE_DRIVE_SCOPE=drive -export RCLONE_CONFIG_GDRIVE_TYPE=drive - -echo ${GOOGLE_SERVICE_JSON} | jq . > ${RCLONE_DRIVE_SERVICE_ACCOUNT_FILE} - -# Little check to make sure we can list from google drive -rclone ls gdrive: - -TEST_DATA_FOLDER=${GITHUB_WORKSPACE}/echopype/test_data -if [ -d $TEST_DATA_FOLDER ] -then - echo "Removing old test data" - rm -rf $TEST_DATA_FOLDER - echo "Copying new test data from google drive" - rclone copy gdrive: $TEST_DATA_FOLDER - echo "Done" - - chmod -R ugoa+w $TEST_DATA_FOLDER - ls -lah $TEST_DATA_FOLDER -else - echo "${TEST_DATA_FOLDER} not found" -fi diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index 64bf8e36c..6912fdb84 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -3,13 +3,16 @@ name: build on: push: branches: - - dev - main paths-ignore: ["**/docker.yaml"] workflow_dispatch: env: NUM_WORKERS: 2 + USE_POOCH: "True" + ECHOPYPE_DATA_VERSION: v0.11.1a2 + ECHOPYPE_DATA_BASEURL: https://github.com/OSOceanAcoustics/echopype/releases/download/{version}/ + XDG_CACHE_HOME: ${{ github.workspace }}/.cache jobs: test: @@ -20,9 +23,10 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.10", "3.11", "3.12"] + python-version: ["3.11", "3.12", "3.13"] runs-on: [ubuntu-latest] experimental: [false] + services: # TODO: figure out how to update tag when there's a new one minio: @@ -30,42 +34,122 @@ jobs: ports: - 9000:9000 httpserver: - image: cormorack/http:latest + image: httpd:2.4 ports: - 8080:80 + steps: - name: Checkout repo - uses: actions/checkout@v4 + uses: actions/checkout@v6 with: fetch-depth: 0 # Fetch all history for all branches and tags. - - name: Set environment variables + + - name: Free disk space + if: runner.os == 'Linux' run: | - echo "PYTHON_VERSION=${{ matrix.python-version }}" >> $GITHUB_ENV + set -eux + echo "Before cleanup:" + df -h + for p in /usr/share/dotnet /opt/ghc /usr/local/lib/android /opt/hostedtoolcache/CodeQL; do + if [ -d "$p" ]; then + echo "Removing $p (size: $(du -sh "$p" | cut -f1))" + sudo rm -rf "$p" + else + echo "Not present: $p" + fi + done + sudo apt-get clean + sudo rm -rf /var/lib/apt/lists/* + echo "After cleanup:" + df -h + - name: Set up Python - uses: actions/setup-python@v5.5.0 + uses: actions/setup-python@v6.1.0 with: python-version: ${{ matrix.python-version }} + + - name: System usage (baseline) + if: runner.os == 'Linux' + run: | + set -eux + echo "== Disk ==" + df -h + echo "== Memory ==" + free -h + echo "== Swap ==" + swapon --show || true + echo "== Top memory processes ==" + ps -eo pid,ppid,cmd,%mem,%cpu --sort=-%mem | head -n 20 || true + - name: Upgrade pip run: python -m pip install --upgrade pip + + - name: Set environment variables + run: | + echo "PYTHON_VERSION=${{ matrix.python-version }}" >> $GITHUB_ENV + - name: Remove docker-compose python run: sed -i "/docker-compose/d" requirements-dev.txt + - name: Install dev tools run: python -m pip install -r requirements-dev.txt + + - name: Cache echopype test data + uses: actions/cache@v5 + with: + path: ${{ env.XDG_CACHE_HOME }}/echopype + key: ep-data-${{ env.ECHOPYPE_DATA_VERSION }} + + - name: Ensure pooch (for GitHub assets) + run: python -m pip install "pooch>=1.8" + + - name: Install error reporter + if: ${{ matrix.python-version == '3.13' }} + run: python -m pip install pytest-github-actions-annotate-failures + - name: Install echopype run: python -m pip install -e ".[plot]" + - name: Print installed packages run: python -m pip list - - name: Copying test data to services - shell: bash -l {0} + + - name: System usage (after install) + if: runner.os == 'Linux' + run: | + set -eux + df -h + free -h + du -h -d 2 "${{ env.XDG_CACHE_HOME }}" || true + ps -eo pid,ppid,cmd,%mem,%cpu --sort=-%mem | head -n 20 || true + + - name: Pre-fetch Pooch test data + run: | + pytest --collect-only -q + + - name: Copy test data to MinIO and HTTP server run: | python .ci_helpers/docker/setup-services.py --deploy --data-only --http-server ${{ job.services.httpserver.id }} # Check data endpoint curl http://localhost:8080/data/ - - name: Running All Tests + + - name: System usage (after services) + if: runner.os == 'Linux' + run: | + set -eux + df -h + free -h + du -h -d 2 "${{ env.XDG_CACHE_HOME }}" || true + docker system df || true + ps -eo pid,ppid,cmd,%mem,%cpu --sort=-%mem | head -n 20 || true + + + - name: Running all tests shell: bash -l {0} run: | - pytest -vv -rx --numprocesses=${{ env.NUM_WORKERS }} --max-worker-restart=3 --cov=echopype --cov-report=xml --log-cli-level=WARNING --disable-warnings + pytest -vvv -rx --numprocesses=${{ env.NUM_WORKERS }} --max-worker-restart=3 \ + --cov=echopype --cov-report=xml --log-cli-level=WARNING --disable-warnings + - name: Upload code coverage to Codecov uses: codecov/codecov-action@v5 with: @@ -74,3 +158,18 @@ jobs: env_vars: RUNNER_OS,PYTHON_VERSION name: codecov-umbrella fail_ci_if_error: false + + - name: System usage (final) + if: always() && runner.os == 'Linux' + run: | + set -eux + echo "== Disk ==" + df -h + echo "== Memory ==" + free -h + echo "== Cache dir ==" + du -h -d 2 "${{ env.XDG_CACHE_HOME }}" || true + echo "== Docker ==" + docker system df || true + echo "== Top memory processes ==" + ps -eo pid,ppid,cmd,%mem,%cpu --sort=-%mem | head -n 30 || true diff --git a/.github/workflows/docker.yaml b/.github/workflows/docker.yaml index 55ee4b53e..5b5ade877 100644 --- a/.github/workflows/docker.yaml +++ b/.github/workflows/docker.yaml @@ -15,31 +15,53 @@ jobs: runs-on: ubuntu-22.04 strategy: matrix: - image_name: ["minioci", "http"] + image_name: ["minioci"] steps: - name: Checkout - uses: actions/checkout@v4 + uses: actions/checkout@v6 + + - run: df -h + - name: Free disk space + uses: jlumbroso/free-disk-space@main + with: + android: true + dotnet: true + haskell: true + tool-cache: true + swap-storage: true + - run: df -h + + - name: Docker prune + run: | + docker buildx prune -af --filter until=24h || true # keep super-recent cache + docker system prune -af || true + - name: Retrieve test data if: matrix.image_name == 'http' uses: ./.github/actions/gdrive-rclone env: GOOGLE_SERVICE_JSON: ${{ secrets.GOOGLE_SERVICE_JSON }} ROOT_FOLDER_ID: ${{ secrets.TEST_DATA_FOLDER_ID }} + - name: Set Docker Image Spec run: | DATE_TAG="$( date -u '+%Y.%m.%d' )" IMAGE_SPEC="${{ env.DOCKER_ORG }}/${{ matrix.image_name }}" echo "IMAGE_SPEC=${IMAGE_SPEC}" >> $GITHUB_ENV echo "DATE_TAG=${DATE_TAG}" >> $GITHUB_ENV + - name: Set up QEMU uses: docker/setup-qemu-action@v3 + - name: Set up Docker Buildx uses: docker/setup-buildx-action@v3 + - name: Login to DockerHub uses: docker/login-action@v3 with: username: ${{ secrets.DOCKER_USERNAME }} password: ${{ secrets.DOCKER_TOKEN }} + - name: Build and push id: docker_build_push uses: docker/build-push-action@v6 diff --git a/.github/workflows/ep-install.yaml b/.github/workflows/ep-install.yaml index 9d5964ecb..d4139574b 100644 --- a/.github/workflows/ep-install.yaml +++ b/.github/workflows/ep-install.yaml @@ -18,7 +18,7 @@ jobs: os: ["ubuntu-latest", "macos-latest", "windows-latest"] steps: - name: Checkout repo - uses: actions/checkout@v4 + uses: actions/checkout@v6 - name: Setup miniconda uses: conda-incubator/setup-miniconda@v3 with: diff --git a/.github/workflows/packit.yaml b/.github/workflows/packit.yaml index b439dab56..1f0f12237 100644 --- a/.github/workflows/packit.yaml +++ b/.github/workflows/packit.yaml @@ -14,13 +14,13 @@ jobs: if: github.repository == 'OSOceanAcoustics/echopype' steps: - name: Checkout - uses: actions/checkout@v4 + uses: actions/checkout@v6 with: # fetch all history so that setuptools-scm works fetch-depth: 0 - name: Set up Python - uses: actions/setup-python@v5.5.0 + uses: actions/setup-python@v6.1.0 with: python-version: 3.12 @@ -42,7 +42,7 @@ jobs: echo "" echo "Generated files:" ls -lh dist/ - - uses: actions/upload-artifact@v4 + - uses: actions/upload-artifact@v6 with: name: releases path: dist @@ -52,11 +52,11 @@ jobs: needs: build-artifact runs-on: ubuntu-22.04 steps: - - uses: actions/setup-python@v5.5.0 + - uses: actions/setup-python@v6.1.0 name: Install Python with: python-version: 3.12 - - uses: actions/download-artifact@v4 + - uses: actions/download-artifact@v7 with: name: releases path: dist diff --git a/.github/workflows/pr.yaml b/.github/workflows/pr.yaml index 4af755a47..46cdec27e 100644 --- a/.github/workflows/pr.yaml +++ b/.github/workflows/pr.yaml @@ -3,9 +3,18 @@ name: Test PR on: pull_request: paths-ignore: ["**/docker.yaml", "docs"] + push: + branches: + - ci/* + workflow_dispatch: + env: NUM_WORKERS: 2 + USE_POOCH: "True" + ECHOPYPE_DATA_VERSION: v0.11.1a2 + ECHOPYPE_DATA_BASEURL: https://github.com/OSOceanAcoustics/echopype/releases/download/{version}/ + XDG_CACHE_HOME: ${{ github.workspace }}/.cache jobs: test: @@ -16,9 +25,11 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.10", "3.11", "3.12"] + python-version: ["3.11", "3.12", "3.13"] runs-on: [ubuntu-latest] experimental: [false] + + services: # TODO: figure out how to update tag when there's a new one minio: @@ -26,57 +37,112 @@ jobs: ports: - 9000:9000 httpserver: - image: cormorack/http:latest + image: httpd:2.4 ports: - 8080:80 + steps: - name: Checkout repo - uses: actions/checkout@v4 + uses: actions/checkout@v6 with: fetch-depth: 0 # Fetch all history for all branches and tags. + + - name: Free disk space + if: runner.os == 'Linux' + run: | + set -eux + echo "Before cleanup:" + df -h + for p in /usr/share/dotnet /opt/ghc /usr/local/lib/android /opt/hostedtoolcache/CodeQL; do + if [ -d "$p" ]; then + echo "Removing $p (size: $(du -sh "$p" | cut -f1))" + sudo rm -rf "$p" + else + echo "Not present: $p" + fi + done + sudo apt-get clean + sudo rm -rf /var/lib/apt/lists/* + echo "After cleanup:" + df -h + + - name: Set up Python - uses: actions/setup-python@v5.5.0 + uses: actions/setup-python@v6.1.0 with: python-version: ${{ matrix.python-version }} + + - name: Upgrade pip run: python -m pip install --upgrade pip + + - name: Set environment variables - run: | - echo "PYTHON_VERSION=${{ matrix.python-version }}" >> $GITHUB_ENV + run: echo "PYTHON_VERSION=${{ matrix.python-version }}" >> $GITHUB_ENV + - name: Remove docker-compose python run: sed -i "/docker-compose/d" requirements-dev.txt + - name: Install dev tools run: python -m pip install -r requirements-dev.txt - # We only want to install this on one run, because otherwise we'll have - # duplicate annotations. + + - name: Cache echopype test data + uses: actions/cache@v5 + with: + path: ${{ env.XDG_CACHE_HOME }}/echopype + key: ep-data-${{ env.ECHOPYPE_DATA_VERSION }} + + - name: Ensure pooch (for GitHub assets) + run: python -m pip install "pooch>=1.8" + - name: Install error reporter - if: ${{ matrix.python-version == '3.12' }} + if: ${{ matrix.python-version == '3.13' }} run: python -m pip install pytest-github-actions-annotate-failures + - name: Install echopype run: python -m pip install -e ".[plot]" + - name: Print installed packages run: python -m pip list - - name: Copying test data to services + + - name: Finding changed files + id: files + uses: tj-actions/changed-files@v47 + + - name: Print Changed files + run: echo "${{ steps.files.outputs.all_changed_files }}" + + - name: Disk usage (before test data) + if: runner.os == 'Linux' + run: | + set -eux + df -h + du -h -d 3 ${{ env.XDG_CACHE_HOME }} || true + + - name: Pre-fetch Pooch test data + run: | + pytest --collect-only -q + + - name: Copy test data to MinIO and HTTP server run: | python .ci_helpers/docker/setup-services.py --deploy --data-only --http-server ${{ job.services.httpserver.id }} # Check data endpoint curl http://localhost:8080/data/ - - name: Finding changed files - id: files - uses: Ana06/get-changed-files@v2.3.0 - with: - format: 'csv' - - name: Print Changed files - run: echo "${{ steps.files.outputs.added_modified_renamed }}" - - name: Running all Tests - if: contains(github.event.pull_request.title, '[all tests ci]') + + - name: Disk usage (after services deploy) + if: runner.os == 'Linux' run: | - pytest -vvv -rx --numprocesses=${{ env.NUM_WORKERS }} --max-worker-restart=3 --cov=echopype --cov-report=xml --log-cli-level=WARNING --disable-warnings - - name: Running Tests - if: ${{ !contains(github.event.pull_request.title, '[all tests ci]') }} + set -eux + df -h + du -h -d 3 ${{ env.XDG_CACHE_HOME }} || true + docker system df || true + + - name: Running all tests run: | - python .ci_helpers/run-test.py --pytest-args="--log-cli-level=WARNING,-vvv,-rx,--numprocesses=${{ env.NUM_WORKERS }},--max-worker-restart=3,--disable-warnings" --include-cov ${{ steps.files.outputs.added_modified_renamed }} + pytest -vvv -rx --numprocesses=${{ env.NUM_WORKERS }} --max-worker-restart=3 \ + --cov=echopype --cov-report=xml --log-cli-level=WARNING --disable-warnings + - name: Upload code coverage to Codecov uses: codecov/codecov-action@v5 with: @@ -85,3 +151,65 @@ jobs: env_vars: RUNNER_OS,PYTHON_VERSION name: codecov-umbrella fail_ci_if_error: false + + test-windows: + name: ${{ matrix.python-version }}--windows + runs-on: windows-latest + strategy: + fail-fast: false + matrix: + python-version: ["3.11", "3.12", "3.13"] + + steps: + - uses: actions/checkout@v6 + with: + fetch-depth: 0 + + - uses: actions/setup-python@v6.1.0 + with: + python-version: ${{ matrix.python-version }} + + - name: Upgrade pip + run: python -m pip install --upgrade pip + + - name: Install deps (dev + pooch + s3) + run: | + python -m pip install -r requirements-dev.txt + python -m pip install "pooch>=1.8" fsspec s3fs + python -m pip install pytest pytest-cov pytest-xdist pytest-mock + + - name: Install error reporter + if: ${{ matrix.python-version == '3.13' }} + run: python -m pip install pytest-github-actions-annotate-failures + + - name: Install echopype + run: python -m pip install -e ".[plot]" + + - name: Pre-fetch all Pooch assets + run: python -m pytest --collect-only -q + + - name: Start local services + run: python .ci_helpers/setup-services-windows.py start + + - name: Running all Tests + shell: pwsh + env: + PY_COLORS: "1" + run: > + python -m pytest -vvv -rx + --numprocesses $env:NUM_WORKERS --max-worker-restart 3 + --cov echopype --cov-report xml + --log-cli-level WARNING --disable-warnings + + - name: Upload code coverage + uses: codecov/codecov-action@v5 + with: + file: ./coverage.xml + flags: unittests + env_vars: RUNNER_OS + name: codecov-windows + fail_ci_if_error: false + + - name: Teardown services + if: always() + run: python .ci_helpers/setup-services-windows.py stop diff --git a/.github/workflows/pypi.yaml b/.github/workflows/pypi.yaml index ab9ac6c47..7c4d20319 100644 --- a/.github/workflows/pypi.yaml +++ b/.github/workflows/pypi.yaml @@ -18,13 +18,13 @@ jobs: if: github.repository == 'OSOceanAcoustics/echopype' || github.event_name == 'workflow_dispatch' steps: - name: Checkout - uses: actions/checkout@v4 + uses: actions/checkout@v6 with: # fetch all history so that setuptools-scm works fetch-depth: 0 - name: Set up Python - uses: actions/setup-python@v5.5.0 + uses: actions/setup-python@v6.1.0 with: python-version: 3.12 @@ -46,7 +46,7 @@ jobs: echo "" echo "Generated files:" ls -lh dist/ - - uses: actions/upload-artifact@v4 + - uses: actions/upload-artifact@v6 with: name: releases path: dist @@ -56,11 +56,11 @@ jobs: needs: build-artifact runs-on: ubuntu-22.04 steps: - - uses: actions/setup-python@v5.5.0 + - uses: actions/setup-python@v6.1.0 name: Install Python with: python-version: 3.12 - - uses: actions/download-artifact@v4 + - uses: actions/download-artifact@v7 with: name: releases path: dist @@ -69,7 +69,7 @@ jobs: ls -ltrh ls -ltrh dist - name: Publish to Test PyPI - uses: pypa/gh-action-pypi-publish@v1.12.4 + uses: pypa/gh-action-pypi-publish@v1.13.0 with: password: ${{ secrets.TEST_PYPI_API_TOKEN }} repository_url: https://test.pypi.org/legacy/ @@ -99,11 +99,11 @@ jobs: if: startsWith(github.ref, 'refs/tags') runs-on: ubuntu-latest steps: - - uses: actions/download-artifact@v4 + - uses: actions/download-artifact@v7 with: name: releases path: dist - name: Publish to PyPI - uses: pypa/gh-action-pypi-publish@v1.12.4 + uses: pypa/gh-action-pypi-publish@v1.13.0 with: password: ${{ secrets.PYPI_API_TOKEN }} diff --git a/.github/workflows/windows-utils.yaml b/.github/workflows/windows-utils.yaml index 77d92e788..6c53440ea 100644 --- a/.github/workflows/windows-utils.yaml +++ b/.github/workflows/windows-utils.yaml @@ -28,7 +28,7 @@ jobs: shell: powershell steps: - name: Checkout repo - uses: actions/checkout@v4 + uses: actions/checkout@v6 with: fetch-depth: 0 # Fetch all history for all branches and tags. - name: Set environment variables diff --git a/.github/workflows/windows.yaml b/.github/workflows/windows.yaml index 582ef5bfe..889fb6188 100644 --- a/.github/workflows/windows.yaml +++ b/.github/workflows/windows.yaml @@ -37,7 +37,7 @@ jobs: - 8080:80 steps: - name: Checkout repo - uses: actions/checkout@v4 + uses: actions/checkout@v6 - name: Copying test data to http server run: | rm .\echopype\test_data -r -fo @@ -46,12 +46,12 @@ jobs: # Check data endpoint curl http://localhost:8080/data/ - name: Setup Python - uses: actions/setup-python@v5.5.0 + uses: actions/setup-python@v6.1.0 with: python-version: ${{ matrix.python-version }} architecture: x64 - name: Cache conda - uses: actions/cache@v4.2.3 + uses: actions/cache@v5 env: # Increase this value to reset cache if '.ci_helpers/py{0}.yaml' has not changed CACHE_NUMBER: 0 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b22e750b4..f02d69042 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -7,7 +7,7 @@ exclude: | ) repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v5.0.0 + rev: v6.0.0 hooks: - id: trailing-whitespace - id: end-of-file-fixer @@ -18,13 +18,13 @@ repos: args: ["--autofix", "--indent=2", "--no-sort-keys"] - repo: https://github.com/PyCQA/isort - rev: 6.0.1 + rev: 7.0.0 hooks: - id: isort args: ["--profile", "black", "--filter-files"] - - repo: https://github.com/psf/black - rev: 25.1.0 + - repo: https://github.com/psf/black-pre-commit-mirror + rev: 25.12.0 hooks: - id: black @@ -39,4 +39,4 @@ repos: - id: codespell # Checks spelling in `docs/source` and `echopype` dirs ONLY # Ignores `.ipynb` files and `_build` folders - args: ["-L", "soop,anc,indx", "--skip=*.ipynb,docs/source/_build,echopype/test_data", "-w", "docs/source", "echopype"] + args: ["-L", "soop,anc,indx,weill", "--skip=*.ipynb,docs/source/_build,echopype/test_data", "-w", "docs/source", "echopype"] diff --git a/docs/requirements.txt b/docs/requirements.txt index a91018353..62464c1a6 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -2,6 +2,6 @@ sphinx_rtd_theme sphinx-automodapi sphinx-panels sphinxcontrib-mermaid -jupyter-book +jupyter-book<2.0 numpydoc docutils<0.18 diff --git a/echopype/calibrate/ek80_complex.py b/echopype/calibrate/ek80_complex.py index 417d6812d..01650b5d7 100644 --- a/echopype/calibrate/ek80_complex.py +++ b/echopype/calibrate/ek80_complex.py @@ -26,7 +26,7 @@ def tapered_chirp( f0 = transmit_frequency_start f1 = transmit_frequency_stop - nsamples = int(np.floor(tau * fs)[0]) + nsamples = int(np.floor(tau * np.float32(fs))[0]) t = np.linspace(0, nsamples - 1, num=nsamples) * 1 / fs a = np.pi * (f1 - f0) / tau b = 2 * np.pi * f0 diff --git a/echopype/clean/__init__.py b/echopype/clean/__init__.py index 9ce5d2908..8a7e7182f 100644 --- a/echopype/clean/__init__.py +++ b/echopype/clean/__init__.py @@ -1,4 +1,5 @@ from .api import ( + detect_transient, estimate_background_noise, mask_attenuated_signal, mask_impulse_noise, @@ -12,4 +13,5 @@ "mask_impulse_noise", "mask_transient_noise", "remove_background_noise", + "detect_transient", ] diff --git a/echopype/clean/api.py b/echopype/clean/api.py index 14924a0e3..a0f451619 100644 --- a/echopype/clean/api.py +++ b/echopype/clean/api.py @@ -11,6 +11,8 @@ from ..utils.compute import _lin2log, _log2lin from ..utils.log import _init_logger from ..utils.prov import add_processing_level, echopype_prov_attrs, insert_input_processing_level +from .transient_noise.transient_fielding import transient_noise_fielding +from .transient_noise.transient_matecho import transient_noise_matecho from .utils import ( add_remove_background_noise_attrs, downsample_upsample_along_depth, @@ -507,3 +509,147 @@ def remove_background_noise( ds_Sv = insert_input_processing_level(ds_Sv, input_ds=ds_Sv) return ds_Sv + + +# Registry of supported methods +METHODS_TRANSIENT = { + "fielding": transient_noise_fielding, + "matecho": transient_noise_matecho, +} + + +def detect_transient(ds: xr.Dataset, method: str, params: dict) -> xr.DataArray: + """ + Dispatch transient-noise detection to a chosen method and return a boolean mask. + + This dispatcher forwards ``ds`` and ``params`` to the selected implementation + (e.g. ``"fielding"``, ``"matecho"``). Any optional arguments omitted in + ``params`` are filled by that method’s own defaults. + + Common expectations + ------------------- + - ``ds`` must contain Sv in dB (default variable name ``"Sv"``). + - Sv should include at least ``ping_time`` and a vertical/range coordinate + (e.g. ``range_sample`` or a ``depth``). + - Returned mask is aligned to Sv (same dims/order) with semantics + **True = VALID (keep)**, **False = transient noise**. + + Parameters + ---------- + ds : xr.Dataset + Acoustic dataset with Sv and required coordinates. + method : str + Name of the detection method. Supported: + - ``"fielding"`` → modified Fielding-style deep transient detector + - ``"matecho"`` → Matecho-style column detector using local percentile + - ``"ryan"`` → to be implemented **TODO** + params : dict + Method-specific keyword arguments (see below). Omitted keys fall back to + the method’s defaults. + + Method-specific arguments + ------------------------- + fielding + var_name : str, default "Sv" + Name of the Sv variable (dB). + range_var : str, default "depth" + Name of the vertical coordinate. Can be 1-D (range) or 2-D (time×range). + r0, r1 : float, default 900, 1000 + Upper/lower bounds of the vertical window (m). + n : int, default 30 + Half-width of the temporal neighbourhood (pings) for the block median. + thr : (float, float), default (3, 1) + Two-stage dB thresholds for detection/propagation. + roff : float, default 20 + Stop depth (m) for upward propagation (don’t propagate above this). + jumps : float, default 5 + Vertical step (m) when iteratively moving the window upward. + maxts : float, default -35 + Max 75th-percentile Sv (dB) for a ping to be considered “quiet”. + start : int, default 0 + Number of initial pings treated as uncomputable in the auxiliary mask + (note: they are still kept/VALID unless you decide otherwise). + + matecho + var_name : str, default "Sv" + Name of the Sv variable (dB). + range_var : str, default "depth" + Vertical coordinate; reduced to 1-D per leading dim if 2-D. + time_var : str, default "ping_time" + Time/ping dimension name. + bottom_var : str | None, default None + Optional bottom depth per ping (m). If None/NaN, uses max range. + start_depth : float, default 220 + Top of the vertical detection window (m). + window_meter : float, default 450 + Vertical window thickness (m). + window_ping : int, default 100 + Temporal window width (pings) for local reference stats. + percentile : float, default 25 + Reference percentile (dB) of the local window. + delta_db : float, default 12 + Excess over percentile (dB) required to flag a ping column. + extend_ping : int, default 0 + Optional horizontal dilation (pings) of flagged columns. + min_window : float, default 20 + Minimum usable window height (m); skip if smaller. + + ryan (**TODO**) + + + Returns + ------- + xr.DataArray + Boolean mask aligned to ``ds[var_name]`` (same dims and order), where + **True = VALID (keep)** and **False = transient noise**. The name/attrs + are set by the called method. + + Raises + ------ + ValueError + If ``method`` is not supported. + + Examples + -------- + >>> from echopype.clean import detect_transient + >>> mask_fielding = detect_transient( + ... ds=ds_Sv_trimmed, + ... method="fielding", + ... params={ + ... "var_name": "Sv", + ... "range_var": "depth", + ... "r0": 900, + ... "r1": 1000, + ... "n": 10, + ... "thr": (3, 1), + ... "roff": 20, + ... "jumps": 5, + ... "maxts": -35, + ... "start": 0, + ... }, + ... ) + + or + + >>> from echopype.clean import detect_transient + >>> mask_matecho = detect_transient( + ... ds=ds_Sv, + ... method="matecho", + ... params={ + ... "var_name": "Sv", + ... "range_var": "depth", + ... "bottom_var": None, + ... "start_depth": 700, + ... "window_meter": 300, + ... "window_ping": 50, + ... "percentile": 25, + ... "delta_db": 8, + ... "extend_ping": 0, + ... "min_window": 5, + ... }, + ... ) + """ + if method not in METHODS_TRANSIENT: + raise ValueError(f"Unsupported transient noise removal method: {method}") + + return METHODS_TRANSIENT[method](ds, **params) diff --git a/echopype/clean/transient_noise/transient_fielding.py b/echopype/clean/transient_noise/transient_fielding.py new file mode 100644 index 000000000..cafcfc547 --- /dev/null +++ b/echopype/clean/transient_noise/transient_fielding.py @@ -0,0 +1,235 @@ +import numpy as np +import xarray as xr + +from echopype.utils.compute import _lin2log, _log2lin + + +def _fielding_core_numpy( + Sv_pr, + r, + r0, + r1, + n, + thr, + roff, + jumps=5, + maxts=-35, + start=0, +): + """ + Core Fielding detector (NumPy). + + Parameters + ---------- + Sv_pr : array-like, shape (ping, range) or (range, ping) + Sv in dB. Internally converted to shape (range, ping). + r : array-like, shape (range,) + Vertical coordinate in meters (monotonic). Used to compute steps/indices. + r0, r1 : float + Vertical window bounds (m). If invalid/outside data, returns all-False mask. + n : int + Half-width of temporal neighbourhood (pings). + thr : tuple[float, float] + Thresholds (dB) for decision stages. + roff : float + Stop depth (m) for upward propagation. + jumps : float + Vertical step (m) when moving the window upward. + maxts : float + Max 75th percentile (dB) to treat a ping as “quiet”. + start : int + Number of initial pings flagged uncomputable in aux mask. + + Returns + ------- + mask_bad_full : np.ndarray of bool, shape (ping, range) + True = BAD (transient noise). + mask_aux_full : np.ndarray of bool, shape (ping, range) + True = uncomputable ping (e.g., insufficient context). Not used by the wrapper. + """ + # transpose to (range, ping) to match your original code + Sv = np.asarray(Sv_pr).T # (range, ping) + r = np.asarray(r) + + if r0 > r1: + # same behavior as original: nothing masked if window invalid + mask = np.zeros_like(Sv, dtype=bool) + mask_ = np.zeros_like(Sv, dtype=bool) + return mask.T, mask_.T + + if (r0 > r[-1]) or (r1 < r[0]): + mask = np.zeros_like(Sv, dtype=bool) + mask_ = np.zeros_like(Sv, dtype=bool) + return mask.T, mask_.T + + up = np.argmin(abs(r - r0)) + lw = np.argmin(abs(r - r1)) + rmin = np.argmin(abs(r - roff)) + + dr = float(np.nanmedian(np.diff(r))) + sf = max(1, int(round(jumps / dr))) + + mask = np.zeros_like(Sv, dtype=bool) # True = BAD + mask_ = np.zeros_like(Sv, dtype=bool) # True = "uncomputable" ping + + n_pings = Sv.shape[1] + for j in range(start, n_pings): + if (j - n < 0) or (j + n > n_pings - 1) or np.all(np.isnan(Sv[up:lw, j])): + mask_[:, j] = True + else: + pingmedian = _lin2log(np.nanmedian(_log2lin(Sv[up:lw, j]))) + pingp75 = _lin2log(np.nanpercentile(_log2lin(Sv[up:lw, j]), 75)) + blockmedian = _lin2log(np.nanmedian(_log2lin(Sv[up:lw, j - n : j + n]))) + + if (pingp75 < maxts) and ((pingmedian - blockmedian) > thr[0]): + r0_, r1_ = up - sf, up + while r0_ > rmin: + pingmedian = _lin2log(np.nanmedian(_log2lin(Sv[r0_:r1_, j]))) + blockmedian = _lin2log(np.nanmedian(_log2lin(Sv[r0_:r1_, j - n : j + n]))) + r0_, r1_ = r0_ - sf, r1_ - sf + if (pingmedian - blockmedian) < thr[1]: + break + mask[r0_:, j] = True + + # restore to (ping, range) for xarray core-dims ("ping_time","range_sample") + # also: pad the first `start` pings with False (good) / True (aux) so shape == input + if start > 0: + pad_bad = np.zeros((start, Sv.shape[0]), dtype=bool) # False = keep + pad_aux = np.ones((start, Sv.shape[0]), dtype=bool) # True = uncomputable + mask_bad_full = np.vstack([pad_bad, mask.T[:, : Sv.shape[1] - start]]) + mask_aux_full = np.vstack([pad_aux, mask_.T[:, : Sv.shape[1] - start]]) + else: + mask_bad_full = mask.T + mask_aux_full = mask_.T + + return mask_bad_full, mask_aux_full + + +def transient_noise_fielding( + ds_Sv: xr.Dataset, + var_name: str = "Sv", + range_var: str = "depth", + r0: float = 900, + r1: float = 1000, + n: int = 30, + thr=(3, 1), + roff: float = 20, + jumps: float = 5, + maxts: float = -35, + start: int = 0, +) -> xr.DataArray: + """ + Transient noise detector modified from the "fielding" + function in `mask_transient.py`, originally written by + Alejandro ARIZA for the Echopy library © 2020. + + Overview + --------- + This algorithm identifies transient noise in echosounder data at deeper + part of the echogram by comparing the echo level of each ping within + its local temporal neighbourhood in a deep water window. + It operates in linear Sv space and in a two-stage process: + + 1) Depth window: In the specified depth interval (`r0`-`r1`, e.g., 900–1000 m), + compute the ping median and the median over neighbouring pings. + If the ping’s deep-window 75th percentile is below `maxts` (i.e., + the window is not broadly high), and the ping median exceeds the + neighbourhood median by more than `thr[0]`, mark the ping as potentially + transient. + + 2) Upward propagation: Move the vertical window upward in fixed steps + (`jumps`, e.g., 5 m). Continue masking shallower ranges until the difference + between the ping median and the median over neighbouring pings drops + below `thr[1]`. This limits the mask to only the part of the water column + affected by the transient. + + Parameters + ---------- + ds_Sv : xr.Dataset + Dataset containing `var_name` (Sv, in dB) and `range_var`. + var_name : str, default "Sv" + Name of the Sv DataArray. Must have dims (..., "ping_time", "range_sample"). + range_var : str, default "depth" + Name of the vertical coordinate. Can be 1-D over ("range_sample") or 2-D over + ("ping_time","range_sample"); the wrapper reduces it to a 1-D vector per + leading dimension (e.g., channel) by selecting the first ping. + r0, r1 : float + Upper/lower bounds of the vertical window (meters). If the window is invalid + or outside the data range, nothing is masked. + n : int + Half-width of the temporal neighbourhood (pings) used to compute the block median. + thr : tuple[float, float], default (3, 1) + Thresholds (dB) used in the two-stage decision. + roff : float + Minimum depth (m) to which the masking can propagate upward (stop depth). + jumps : float, default 5 + Vertical step (m) used when iteratively moving the window upward. + maxts : float, default -35 + Maximum allowable 75th-percentile Sv (dB) to consider a ping “quiet enough.” + start : int, default 0 + Number of initial pings to mark as uncomputable in the auxiliary mask (internal). + + Returns + ------- + xr.DataArray (bool) + Boolean mask aligned to `ds_Sv[var_name]` with the **same dims and order**, + where **True = VALID (keep)** and **False = transient noise**. + Name: "fielding_mask_valid". Dtype: bool. + + Examples, to be used with dispatcher + -------- + >>> mask_fielding = mask_transient_noise_dispatch( + ds=ds_Sv, + method="fielding", + params={ + "var_name": "Sv", + "range_var": "depth", + "r0": 900, + "r1": 1000, + "n": 10, + "thr": (3, 1), + "roff": 20, + "jumps": 5, + "maxts": -35, + "start": 0, + }, + ) + """ + + if var_name not in ds_Sv: + raise ValueError(f"{var_name!r} not found in Dataset.") + if range_var not in ds_Sv: + raise ValueError(f"{range_var!r} not found in Dataset.") + + Sv_da = ds_Sv[var_name] + + # 1D range vector along 'range_sample' + r_da = ds_Sv[range_var] + if {"ping_time", "range_sample"}.issubset(r_da.dims): + r_1d = r_da.isel(ping_time=0) + elif (r_da.ndim == 1) and ("range_sample" in r_da.dims): + r_1d = r_da + else: + raise ValueError(f"Cannot infer 1D '{range_var}' from dims {r_da.dims}.") + + mask_bad, mask_aux = xr.apply_ufunc( + _fielding_core_numpy, + Sv_da, + r_1d, + input_core_dims=[["ping_time", "range_sample"], ["range_sample"]], + output_core_dims=[["ping_time", "range_sample"], ["ping_time", "range_sample"]], + vectorize=True, + dask="parallelized", + output_dtypes=[bool, bool], + kwargs=dict(r0=r0, r1=r1, n=n, thr=thr, roff=roff, jumps=jumps, maxts=maxts, start=start), + ) + + # Flip to True = VALID for apply_mask + mask_valid = ( + (~mask_bad) + .astype(bool) + .rename("fielding_mask_valid") + .assign_attrs({"meaning": "True = VALID (False = transient noise)"}) + ) + # Ensure dims order matches `Sv_da` + return mask_valid.transpose(*Sv_da.dims) diff --git a/echopype/clean/transient_noise/transient_matecho.py b/echopype/clean/transient_noise/transient_matecho.py new file mode 100644 index 000000000..cd6eec810 --- /dev/null +++ b/echopype/clean/transient_noise/transient_matecho.py @@ -0,0 +1,287 @@ +import numpy as np +import xarray as xr +from scipy.ndimage import binary_dilation + +from echopype.utils.compute import _lin2log, _log2lin + + +def _matecho_core_numpy( + Sv, + r, + bottom_depth=None, + start_depth=220, + window_meter=450, + window_ping=100, + percentile=25, + delta_db=12, + extend_ping=0, + min_window=20, +): + """ + Matecho-style transient detector (NumPy core operating on columns). + + Parameters + ---------- + Sv : np.ndarray, shape (range, ping) or (ping, range) + Sv in dB. Internally converted to (range, ping). + r : array-like, shape (range,) + Vertical coordinate (meters). Must match Sv’s range axis length. + bottom_depth : array-like or None, shape (ping,), optional + Bottom depth per ping (meters). If None/NaN, set to r[-1] per ping. + start_depth : float + Top of vertical detection window (m). + window_meter : float + Vertical window thickness (m). + window_ping : int + Temporal window width (pings) centered on each ping for local stats. + percentile : float + Reference percentile (dB) computed from the local window. + delta_db : float + Excess over percentile (dB) required to flag a ping. + extend_ping : int + Horizontal dilation (pings) for flagged columns. + min_window : float + Minimum usable window height (m). Skip if smaller. + + Returns + ------- + mask_bad : np.ndarray of bool, shape (range, ping) + True = BAD (transient noise) — columns for flagged pings. + aux_2d : np.ndarray of bool, shape (range, ping) + Reserved for diagnostics (currently all False). + """ + + n_ping = Sv.shape[1] + r = np.asarray(r) + depth_mask = (r >= start_depth) & (r <= start_depth + window_meter) + + # Sv must be (range, ping) + if Sv.shape[0] != len(r): + if Sv.shape[1] == len(r): + Sv = Sv.T + else: + raise ValueError(f"Mismatch Sv {Sv.shape} vs range {len(r)}") + + # bottom_depth + if bottom_depth is None: + bottom_depth = np.full(n_ping, r[-1], dtype=float) + else: + bottom_depth = np.array(bottom_depth, dtype=float, copy=True) + bottom_depth[np.isnan(bottom_depth)] = r[-1] + + pings_bad = np.zeros(n_ping, dtype=bool) + + for j in range(n_ping): + j0 = max(0, j - window_ping // 2) + j1 = min(n_ping, j + window_ping // 2) + local_bottom = np.min(bottom_depth[j0:j1]) + + refined_mask = depth_mask & (r < local_bottom) + if not np.any(refined_mask): + continue + + H = (r[1] - r[0]) * np.sum(refined_mask) + if H < min_window: + continue + + sv_window = Sv[refined_mask, j0:j1] + sv_ping = Sv[refined_mask, j] + + sv_window_flat = sv_window[~np.isnan(sv_window)] + if sv_window_flat.size == 0: + continue + + pctl_val = np.percentile(sv_window_flat, percentile) # could compute in linear + sv_ping_lin = _log2lin(sv_ping) + sv_ping_mean_db = _lin2log(np.nanmean(sv_ping_lin)) + + if sv_ping_mean_db > pctl_val + delta_db: + pings_bad[j] = True + + if extend_ping > 0: + structure = np.ones(2 * extend_ping + 1, dtype=bool) + pings_bad = binary_dilation(pings_bad, structure=structure) + + mask_bad = np.zeros_like(Sv, dtype=bool) # (range, ping) + mask_bad[:, pings_bad] = True + aux_2d = np.zeros_like(mask_bad, dtype=bool) + return mask_bad, aux_2d # True = BAD + + +def transient_noise_matecho( + ds: xr.Dataset, + var_name: str = "Sv", + range_var: str = "depth", + time_var: str = "ping_time", + bottom_var: str | None = None, + start_depth: float = 220, + window_meter: float = 450, + window_ping: int = 100, + percentile: float = 25, + delta_db: float = 12, + extend_ping: int = 0, + min_window: float = 20, +) -> xr.DataArray: + """ + Transient noise detector modified from the "DeepSpikeDetection" + function in `DeepSpikeDetection.m`, originally written by + Yannick PERROT for the Matecho open-source tool for processing + fisheries acoustics data © 2018. + + Perrot, Y., Brehmer, P., Habasque, J., Roudaut, G., Behagle, N., Sarré, A. + and Lebourges-Dhaussy, A., 2018. Matecho: an open-source tool for processing + fisheries acoustics data. Acoustics Australia, 46(2), pp.241-248. + + Overview + -------- + Flags entire pings as transient when, within a deep window, the ping’s + mean Sv (computed in linear units, then converted back to dB) exceeds a + local reference percentile by `delta_db`. + + 1) Depth window: Use a vertical slice from `start_depth` to + `start_depth + window_meter`, limited by a local bottom (if provided; + otherwise will use the maximum depth of the echogram). + Skip if usable height < `min_window`. + + 2) Local reference: For ping j, form a temporal neighborhood + [j - window_ping/2, j + window_ping/2] and compute the chosen `percentile` + over that neighborhood within the deep window. + + 3) Mean Sv within the depth window (`ping_mean_db`): + Compute the mean Sv (in the linear domain and converted back to dB). + + 4) Decision: If `ping_mean_db > percentile + delta_db`, mark ping j as BAD. + Optionally dilate flagged pings horizontally by `extend_ping` + (binary dilation). + + This function prepares the vertical coordinate (and optional bottom), then + calls a NumPy core via `xarray.apply_ufunc`, vectorized across leading dims + (e.g., `channel`). + + Parameters + ---------- + ds : xr.Dataset + Dataset containing `var_name` (Sv in dB) and `range_var`. + var_name : str, default "Sv" + Name of the Sv DataArray. Must include dims (..., "ping_time", "range_sample"). + range_var : str, default "depth" + Vertical coordinate. Can be 1-D or 2-D; reduced to a 1-D vector per leading dim. + time_var : str, default "ping_time" + Time/ping dimension name. + bottom_var : str | None, default None + Name of a 1-D bottom-depth-per-ping variable (meters). If None or NaN, the + maximum range is used. **Currently not plumbed through**; future work. + start_depth : float, default 220 + Top of the vertical detection window (m). + window_meter : float, default 450 + Vertical window thickness (m). + window_ping : int, default 100 + Temporal window width (pings) for the local reference statistics. + percentile : float, default 25 + Reference percentile (in dB space) of the local window. + delta_db : float, default 12 + Threshold (dB) added to the percentile for flagging. + extend_ping : int, default 0 + Optional horizontal dilation in pings applied to flagged columns. + min_window : float, default 20 + Minimum usable vertical window height (m); skip if smaller. + + Returns + ------- + xr.DataArray (bool) + Boolean mask aligned to `ds[var_name]` with the **same dims and order**, + where **True = VALID (keep)** and **False = transient noise**. + Name: "matecho_mask_valid". Dtype: bool. + + Examples to use with dispatcher + -------- + >>> mask_matecho = mask_transient_noise_dispatch( + ds=ds_Sv, + method="matecho", + params={ + "var_name": "Sv", + "range_var": "depth", + "time_var": "ping_time", + "bottom_var": None, + "start_depth": 700, + "window_meter": 300, + "window_ping": 50, + "percentile": 25, + "delta_db": 8, + "extend_ping": 0, + "min_window": 5, + }, + ) + """ + + if var_name not in ds: + raise ValueError(f"{var_name!r} not found.") + if range_var not in ds: + raise ValueError(f"{range_var!r} not found.") + if time_var not in ds[var_name].dims: + raise ValueError(f"{time_var!r} must be a dim of {var_name!r}.") + + Sv_da = ds[var_name] + + rng_dim = None + for cand in (range_var, "range_sample"): + if cand in Sv_da.dims: + rng_dim = cand + break + if rng_dim is None: + raise ValueError(f"No range dim in {var_name!r} dims {Sv_da.dims}") + + # per-channel depth vector (channel, range_sample) -> drop time + r_1d = ds[range_var].isel({time_var: 0}) # keeps 'channel' + + # 1D time vector + t_1d = ds[var_name][time_var] # ('ping_time',) + if t_1d.ndim != 1: + t_1d = t_1d.squeeze() + + # if provided, use ds[bottom_var] (1D per ping); else NaN → core falls back to r[-1] + if bottom_var is not None and bottom_var in ds: + bottom_1d = ds[bottom_var].transpose(..., time_var).astype(float) + else: + bottom_1d = xr.DataArray(np.full(ds[var_name][time_var].shape, np.nan), dims=(time_var,)) + + # ensure Sv’s last two dims are (range, time) for the core (keep leading dims like channel) + Sv_core = Sv_da.transpose(..., rng_dim, time_var) + + mask_bad, _ = xr.apply_ufunc( + _matecho_core_numpy, + Sv_core, + r_1d, + bottom_1d, + input_core_dims=[ + [rng_dim, time_var], + [rng_dim], + [time_var], + ], + output_core_dims=[ + [rng_dim, time_var], + [rng_dim, time_var], + ], + vectorize=True, + dask="parallelized", + output_dtypes=[bool, bool], + kwargs=dict( + start_depth=start_depth, + window_meter=window_meter, + window_ping=window_ping, + percentile=percentile, + delta_db=delta_db, + extend_ping=extend_ping, + min_window=min_window, + ), + ) + + mask_valid = ( + (~mask_bad) + .astype(bool) + .rename("matecho_mask_valid") + .assign_attrs({"meaning": "True = VALID (False = transient noise)"}) + ) + + # Return with SAME dims/order as input Sv (e.g., ('channel','ping_time','range_sample')) + return mask_valid.transpose(*Sv_da.dims) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 54c1e2dcf..54e3b71f9 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -137,13 +137,16 @@ def compute_MVBS( **flox_kwargs, ) + # Generalize the first dimension name to support multiple like channel and frequency_nominal + dim_0 = list(raw_MVBS.sizes.keys())[0] + # create MVBS dataset # by transforming the binned dimensions to regular coords ds_MVBS = xr.Dataset( - data_vars={"Sv": (["channel", "ping_time", range_var], raw_MVBS["Sv"].data)}, + data_vars={"Sv": ([dim_0, "ping_time", range_var], raw_MVBS["Sv"].data)}, coords={ "ping_time": np.array([v.left for v in raw_MVBS.ping_time_bins.values]), - "channel": raw_MVBS.channel.values, + dim_0: getattr(raw_MVBS, dim_0).values, range_var: np.array([v.left for v in raw_MVBS[f"{range_var}_bins"].values]), }, ) @@ -206,6 +209,7 @@ def compute_MVBS( prov_dict["processing_function"] = "commongrid.compute_MVBS" ds_MVBS = ds_MVBS.assign_attrs(prov_dict) ds_MVBS["frequency_nominal"] = ds_Sv["frequency_nominal"] # re-attach frequency_nominal + ds_MVBS["channel"] = ds_Sv["channel"] # re-attach channel ds_MVBS = insert_input_processing_level(ds_MVBS, input_ds=ds_Sv) diff --git a/echopype/commongrid/utils.py b/echopype/commongrid/utils.py index 962820db9..d831c39b1 100644 --- a/echopype/commongrid/utils.py +++ b/echopype/commongrid/utils.py @@ -35,6 +35,7 @@ def compute_raw_MVBS( with coordinates ``channel``, ``ping_time``, and ``range_sample`` at bare minimum. Or this can contain ``Sv`` and ``depth`` data with similar coordinates. + ``frequency_nominal`` is supported as an alternative to ``channel`` range_interval: pd.IntervalIndex or np.ndarray 1D array or interval index representing the bins required for ``range_var`` @@ -530,6 +531,8 @@ def _groupby_x_along_channels( For NASC computatioon this must contain ``Sv`` and ``depth`` data with coordinates ``channel``, ``distance_nmi``, and ``range_sample``. + + ``frequency_nominal`` is supported as an alternative to ``channel`` range_interval: pd.IntervalIndex or np.ndarray 1D array or interval index representing the bins required for ``range_var`` @@ -604,12 +607,13 @@ def _groupby_x_along_channels( f"The ```{array_name}``` coordinate array contain NaNs. {aggregation_msg}" ) - # reduce along ping_time or distance_nmi - # and echo_range or depth - # by binning and averaging + # Use the first dimension as the grouping dimension for generality + dim_0 = list(ds_Sv.sizes.keys())[0] + + # bin and average along ping_time or distance_nmi and echo_range or depth sv_mean = xarray_reduce( sv, - ds_Sv["channel"], + ds_Sv[dim_0], # generic: not always 'channel' ds_Sv[x_var], ds_Sv[range_var], expected_groups=(None, x_interval, range_interval), diff --git a/echopype/consolidate/api.py b/echopype/consolidate/api.py index 26f4ea6ba..81089d047 100644 --- a/echopype/consolidate/api.py +++ b/echopype/consolidate/api.py @@ -26,6 +26,7 @@ logger = _init_logger(__name__) POSITION_VARIABLES = ["latitude", "longitude"] +SUPPORTED_DIM_0_NAMES = ["channel", "frequency_nominal"] def swap_dims_channel_frequency(ds: Union[xr.Dataset, str, pathlib.Path]) -> xr.Dataset: @@ -64,6 +65,28 @@ def swap_dims_channel_frequency(ds: Union[xr.Dataset, str, pathlib.Path]) -> xr. ) +def get_dim_0(ds: Union[xr.Dataset, str, pathlib.Path]) -> str: + """ + Get the name of the first dimension of the dataset. + + Parameters + ---------- + ds : Union[xr.Dataset, str, pathlib.Path] + The input dataset. + + Returns + ------- + str + The name of the first dimension. + """ + if list(ds.sizes.keys())[0] in SUPPORTED_DIM_0_NAMES: + return list(ds.sizes.keys())[0] + else: + raise ValueError( + f"The first dimension of the dataset must be one of {SUPPORTED_DIM_0_NAMES}." + ) + + @add_processing_level("L2A") def add_depth( ds: Union[xr.Dataset, str, pathlib.Path], @@ -205,12 +228,12 @@ def add_depth( # Compute echo range scaling in EK systems using platform angle data echo_range_scaling = ek_use_platform_angles(echodata["Platform"], ds["ping_time"]) elif use_beam_angles: + dim_0 = get_dim_0(ds) # Identify beam group name by checking channel values of `ds` - if echodata["Sonar/Beam_group1"]["channel"].equals(ds["channel"]): + if echodata["Sonar/Beam_group1"][dim_0].equals(ds[dim_0]): beam_group_name = "Beam_group1" else: beam_group_name = "Beam_group2" - # Compute echo range scaling in EK systems using beam angle data echo_range_scaling = ek_use_beam_angles(echodata[f"Sonar/{beam_group_name}"]) @@ -456,12 +479,8 @@ def add_splitbeam_angle( # and obtain the echodata group path corresponding to encode_mode ed_beam_group = retrieve_correct_beam_group(echodata, waveform_mode, encode_mode) - # check that source_Sv at least has a channel dimension - if "channel" not in source_Sv.variables: - raise ValueError("The input source_Sv Dataset must have a channel dimension!") - - # Select ds_beam channels from source_Sv - ds_beam = echodata[ed_beam_group].sel(channel=source_Sv["channel"].values) + dim_0 = get_dim_0(source_Sv) + ds_beam = echodata[ed_beam_group].sel({dim_0: source_Sv[dim_0].values}) # Assemble angle param dict angle_param_list = [ @@ -478,10 +497,9 @@ def add_splitbeam_angle( raise ValueError(f"source_Sv does not contain the necessary parameter {p_name}!") # fail if source_Sv and ds_beam do not have the same lengths - # for ping_time, range_sample, and channel + # for dim_0, ping_time, range_sample same_size_lens = [ - ds_beam.sizes[dim] == source_Sv.sizes[dim] - for dim in ["channel", "ping_time", "range_sample"] + ds_beam.sizes[dim] == source_Sv.sizes[dim] for dim in [dim_0, "ping_time", "range_sample"] ] if not same_size_lens: raise ValueError( @@ -501,7 +519,7 @@ def add_splitbeam_angle( if pulse_compression: # with pulse compression # put receiver fs into the same dict for simplicity pc_params = get_filter_coeff( - echodata["Vendor_specific"].sel(channel=source_Sv["channel"].values) + echodata["Vendor_specific"].sel({dim_0: source_Sv[dim_0].values}) ) pc_params["receiver_sampling_frequency"] = source_Sv["receiver_sampling_frequency"] theta, phi = get_angle_complex_samples(ds_beam, angle_params, pc_params) diff --git a/echopype/consolidate/split_beam_angle.py b/echopype/consolidate/split_beam_angle.py index a269d70c4..83160b8c5 100644 --- a/echopype/consolidate/split_beam_angle.py +++ b/echopype/consolidate/split_beam_angle.py @@ -210,18 +210,19 @@ def get_angle_complex_samples( else: # beam_type different for some channels, process each channel separately theta, phi = [], [] - for ch_id in bs["channel"].data: + dim_0 = list(bs.sizes.keys())[0] + for ch_id in bs[dim_0].data: theta_ch, phi_ch = _compute_angle_from_complex( - bs=bs.sel(channel=ch_id), + bs=bs.sel({dim_0: ch_id}), # beam_type is not time-varying - beam_type=(ds_beam["beam_type"].sel(channel=ch_id)), + beam_type=(ds_beam["beam_type"].sel({dim_0: ch_id})), sens=[ - angle_params["angle_sensitivity_alongship"].sel(channel=ch_id), - angle_params["angle_sensitivity_athwartship"].sel(channel=ch_id), + angle_params["angle_sensitivity_alongship"].sel({dim_0: ch_id}), + angle_params["angle_sensitivity_athwartship"].sel({dim_0: ch_id}), ], offset=[ - angle_params["angle_offset_alongship"].sel(channel=ch_id), - angle_params["angle_offset_athwartship"].sel(channel=ch_id), + angle_params["angle_offset_alongship"].sel({dim_0: ch_id}), + angle_params["angle_offset_athwartship"].sel({dim_0: ch_id}), ], ) theta.append(theta_ch) @@ -231,7 +232,7 @@ def get_angle_complex_samples( theta = xr.DataArray( data=theta, coords={ - "channel": bs["channel"], + dim_0: bs[dim_0], "ping_time": bs["ping_time"], "range_sample": bs["range_sample"], }, @@ -239,7 +240,7 @@ def get_angle_complex_samples( phi = xr.DataArray( data=phi, coords={ - "channel": bs["channel"], + dim_0: bs[dim_0], "ping_time": bs["ping_time"], "range_sample": bs["range_sample"], }, diff --git a/echopype/convert/parse_ad2cp.py b/echopype/convert/parse_ad2cp.py index c98190495..3cb887af6 100644 --- a/echopype/convert/parse_ad2cp.py +++ b/echopype/convert/parse_ad2cp.py @@ -319,7 +319,7 @@ def timestamp(self) -> np.datetime64: Calculates and returns the timestamp of the packet """ - year = self.data["year"] + 1900 + year = self.data["year"].astype(np.uint16) + 1900 month = self.data["month"] + 1 day = self.data["day"] hour = self.data["hour"] diff --git a/echopype/convert/parse_base.py b/echopype/convert/parse_base.py index 578e99549..ec89baf5b 100644 --- a/echopype/convert/parse_base.py +++ b/echopype/convert/parse_base.py @@ -1,5 +1,6 @@ import datetime import os +import re import sys from collections import defaultdict from typing import Any, Dict, Literal, Optional, Tuple @@ -25,6 +26,15 @@ logger = _init_logger(__name__) +# --- Windows-safe path component sanitizer (also harmless on POSIX) --- +_INVALID_FS_CHARS = r'[<>:"/\\|?*]' + + +def _sanitize_component(s: str) -> str: + """Replace filesystem-invalid characters in path components.""" + return re.sub(_INVALID_FS_CHARS, "_", str(s)) + + class ParseBase: """Parent class for all convert classes.""" @@ -172,9 +182,10 @@ def rectangularize_data( # Setup temp store zarr_store = create_temp_zarr_store() # Setup zarr store - zarr_root = zarr.group( - store=zarr_store, overwrite=True, synchronizer=zarr.ThreadSynchronizer() - ) + # synchronizer not supported yet in zarr v3 + # TODO: add back synchronizer when supported see #1558 + # group(store=zarr_store.root, overwrite=True, synchronizer=zarr.ThreadSynchronizer()) + zarr_root = zarr.group(store=zarr_store.root, overwrite=True) for raw_type in self.raw_types: data_type_shapes = expanded_data_shapes[raw_type] @@ -198,13 +209,12 @@ def _write_to_temp_zarr( chunks: Tuple[int], ) -> dask.array.Array: if shape == arr.shape: - z_arr = zarr_root.array( + z_arr = zarr_root.create_array( name=path, - data=arr, + data=arr.astype(np.float64), fill_value=np.nan, chunks=chunks, - dtype="f8", - write_empty_chunks=False, + config={"write_empty_chunks": False}, ) else: # Figure out current data region @@ -217,7 +227,7 @@ def _write_to_temp_zarr( chunks=chunks, dtype="f8", fill_value=np.nan, # same as float64 - write_empty_chunks=False, + config={"write_empty_chunks": False}, ) # Fill zarr array with actual data @@ -317,32 +327,31 @@ def _parse_and_pad_datagram( # SWAP -------------------------------------------------------------- if data_shape[0] != len(self.ping_time[ch_id]): - # Let's ensure that the ping time dimension is the same - # as the original data shape when written to zarr array - # since this will become the coordinate dimension - # of the channel data when set in set_groups operation + # Ensure ping-time dimension matches original length data_shape = (len(self.ping_time[ch_id]),) + data_shape[1:] - # Write to temp zarr for swap - # then assign the dask array to the dictionary + + # Build safe path components (handles Windows-illegal chars like "|") + safe_ch = _sanitize_component(ch_id) + + # Write to temp zarr for swap, then assign the dask array to the dict if data_type == "complex": # Save real and imaginary components separately ping_data_dict[data_type][ch_id] = {} - # Go through real and imaginary components for complex_part, arr in padded_arr.items(): + safe_part = _sanitize_component(complex_part) d_arr = self._write_to_temp_zarr( arr, zarr_root, - os.path.join(raw_type, data_type, str(ch_id), complex_part), + os.path.join(raw_type, data_type, safe_ch, safe_part), data_shape, chunks, ) ping_data_dict[data_type][ch_id][complex_part] = d_arr - else: d_arr = self._write_to_temp_zarr( padded_arr, zarr_root, - os.path.join(raw_type, data_type, str(ch_id)), + os.path.join(raw_type, data_type, safe_ch), data_shape, chunks, ) diff --git a/echopype/convert/set_groups_ek60.py b/echopype/convert/set_groups_ek60.py index 21c4aced5..d0a12bac0 100644 --- a/echopype/convert/set_groups_ek60.py +++ b/echopype/convert/set_groups_ek60.py @@ -146,7 +146,7 @@ def set_env(self) -> xr.Dataset: ds_env.append(ds_tmp) # Merge data from all channels - ds = xr.merge(ds_env) + ds = xr.merge(ds_env, compat="no_conflicts", join="outer") return set_time_encodings(ds) @@ -308,7 +308,7 @@ def set_platform(self) -> xr.Dataset: # TODO: for current test data we see all # pitch/roll/heave are the same for all freq channels # consider only saving those from the first channel - ds_plat = xr.merge(ds_plat) + ds_plat = xr.merge(ds_plat, compat="no_conflicts", join="outer") ds_plat["channel"] = ds_plat["channel"].assign_attrs( self._varattrs["beam_coord_default"]["channel"] ) diff --git a/echopype/convert/set_groups_ek80.py b/echopype/convert/set_groups_ek80.py index 9fab0260e..930c5b2e4 100644 --- a/echopype/convert/set_groups_ek80.py +++ b/echopype/convert/set_groups_ek80.py @@ -1065,10 +1065,10 @@ def pulse_form_map(pulse_form): def merge_save(ds_combine: List[xr.Dataset], ds_invariant: xr.Dataset) -> xr.Dataset: """Merge data from all complex or all power/angle channels""" # Combine all channels into one Dataset - ds_combine = xr.concat(ds_combine, dim="channel") + ds_combine = xr.concat(ds_combine, dim="channel", join="outer") ds_combine = xr.merge( - [ds_invariant, ds_combine], combine_attrs="override" + [ds_invariant, ds_combine], combine_attrs="override", join="outer" ) # override keeps the Dataset attributes return set_time_encodings(ds_combine) @@ -1342,7 +1342,7 @@ def set_vendor(self) -> xr.Dataset: "long_name" ] = "ID of channels containing broadband calibration information" ds_cal.append(ds_ch) - ds_cal = xr.merge(ds_cal) + ds_cal = xr.merge(ds_cal, join="outer") if "impedance" in ds_cal: ds_cal = ds_cal.rename_vars({"impedance": "impedance_transducer"}) diff --git a/echopype/convert/utils/ek_date_conversion.py b/echopype/convert/utils/ek_date_conversion.py index 78dbc5747..e8384e78f 100644 --- a/echopype/convert/utils/ek_date_conversion.py +++ b/echopype/convert/utils/ek_date_conversion.py @@ -40,7 +40,7 @@ def nt_to_unix(nt_timestamp_tuple, return_datetime=True): The timestamp is a 64bit count of 100ns intervals since the NT epoch broken into two 32bit longs, least significant first: - >>> dt = nt_to_unix((19496896L, 30196149L)) + >>> dt = nt_to_unix((19496896, 30196149)) >>> match_dt = datetime.datetime(2011, 12, 23, 20, 54, 3, 964000, pytz_utc) >>> assert abs(dt - match_dt) <= dt.resolution """ @@ -63,7 +63,7 @@ def unix_to_nt(unix_timestamp): #Simple conversion >>> dt = datetime.datetime(2011, 12, 23, 20, 54, 3, 964000, pytz_utc) - >>> assert (19496896L, 30196149L) == unix_to_nt(dt) + >>> assert (19496896, 30196149) == unix_to_nt(dt) #Converting back and forth between the two standards: >>> orig_dt = datetime.datetime.now(tz=pytz_utc) diff --git a/echopype/echodata/echodata.py b/echopype/echodata/echodata.py index a45d845db..24ac2473f 100644 --- a/echopype/echodata/echodata.py +++ b/echopype/echodata/echodata.py @@ -11,7 +11,7 @@ import numpy as np import xarray as xr from xarray import DataTree, open_datatree, open_groups -from zarr.errors import GroupNotFoundError, PathNotFoundError +from zarr.errors import GroupNotFoundError if TYPE_CHECKING: from ..core import EngineHint, FileFormatHint, PathHint, SonarModelsHint @@ -94,8 +94,8 @@ def cleanup_swap_files(self): ] if len(zarr_stores) > 0: # Grab the first associated file since there is only one unique file - # Can check using zarr_stores[0].path - fs = zarr_stores[0].fs + # Can check using zarr_stores[0] + fs = zarr_stores[0] from ..utils.io import delete_zarr_store for store in zarr_stores: @@ -341,7 +341,7 @@ def __setitem__(self, __key: Optional[str], __newvalue: Any) -> Optional[xr.Data node.dataset = __newvalue return self.__get_dataset(node) except KeyError: - raise GroupNotFoundError(__key) + raise GroupNotFoundError(self.converted_raw_path, __key) else: raise ValueError("Datatree not found!") @@ -472,9 +472,8 @@ def update_platform( ext_var = mappings_expanded[platform_var]["external_var"] platform_var_attrs = platform[platform_var].attrs.copy() - # Create new (replaced) variable using dataset "update" # With update, dropping the variable first is not needed - platform = platform.update({platform_var: (time_dim, ext_ds[ext_var].data)}) + platform.update({platform_var: (time_dim, ext_ds[ext_var].data)}) # Assign attributes to newly created (replaced) variables var_attrs = platform_var_attrs @@ -538,7 +537,7 @@ def _load_file(self, raw_path: "PathHint"): raw_path, group=value["ep_group"], ) - except (OSError, GroupNotFoundError, PathNotFoundError): + except (OSError, GroupNotFoundError): # Skips group not found errors for EK80 and ADCP ... if group == "top" and hasattr(ds, "keywords"): @@ -633,7 +632,6 @@ def to_zarr( overwrite: bool = False, parallel: bool = False, output_storage_options: Dict[str, str] = {}, - consolidated: bool = True, **kwargs, ): """Save content of EchoData to zarr. @@ -652,9 +650,6 @@ def to_zarr( whether or not to use parallel processing. (Not yet implemented) output_storage_options : dict Additional keywords to pass to the filesystem class. - consolidated : bool - Flag to consolidate zarr metadata. - Defaults to ``True`` **kwargs : dict, optional Extra arguments to `xr.Dataset.to_zarr`: refer to xarray's documentation for a list of all possible arguments. @@ -669,7 +664,6 @@ def to_zarr( overwrite=overwrite, parallel=parallel, output_storage_options=output_storage_options, - consolidated=consolidated, **kwargs, ) diff --git a/echopype/mask/__init__.py b/echopype/mask/__init__.py index 24b6c172f..b45d84b33 100644 --- a/echopype/mask/__init__.py +++ b/echopype/mask/__init__.py @@ -1,3 +1,3 @@ -from .api import apply_mask, frequency_differencing +from .api import apply_mask, detect_seafloor, detect_shoal, frequency_differencing -__all__ = ["frequency_differencing", "apply_mask"] +__all__ = ["frequency_differencing", "apply_mask", "detect_seafloor", "detect_shoal"] diff --git a/echopype/mask/api.py b/echopype/mask/api.py index 65b64792b..48c3d88c0 100644 --- a/echopype/mask/api.py +++ b/echopype/mask/api.py @@ -2,17 +2,25 @@ import operator as op import pathlib import sys -from typing import List, Optional, Union +from typing import Dict, List, Optional, Union import dask import dask.array import numpy as np import xarray as xr +# for seafloor detection +from echopype.mask.seafloor_detection.bottom_basic import bottom_basic +from echopype.mask.seafloor_detection.bottom_blackwell import bottom_blackwell + from ..utils.io import validate_source from ..utils.prov import add_processing_level, echopype_prov_attrs, insert_input_processing_level from .freq_diff import _check_freq_diff_source_Sv, _parse_freq_diff_eq +# for schoals detection +from .shoal_detection.shoal_echoview import shoal_echoview +from .shoal_detection.shoal_weill import shoal_weill + # lookup table with key string operator and value as corresponding Python operator str2ops = { ">": op.gt, @@ -658,3 +666,136 @@ def _create_mask(lhs: np.ndarray, diff: float) -> np.ndarray: da = da.assign_attrs(xr_dataarray_attrs) return da + + +# Registry of supported methods for bottom detection +METHODS_BOTTOM = { + "basic": bottom_basic, + "blackwell": bottom_blackwell, +} + + +def detect_seafloor( + ds: xr.Dataset, + method: str, + params: Dict, +) -> xr.DataArray: + """ + Dispatch seafloor detection to a chosen method and return a 1-D bottom line. + + This function forwards ``ds`` and ``params`` to the selected implementation + (e.g., ``"basic"``, ``"blackwell"``). Any optional arguments omitted in + ``params`` are filled by that method's defaults. + + Parameters + ---------- + ds : xr.Dataset + Dataset containing calibrated Sv and required coordinates (at minimum + ``ping_time`` and a vertical coordinate such as ``depth``), plus any + extra variables required by the chosen method. + method : str + Name of the detection method to use. Supported: + - ``"basic"`` → threshold-only detector + - ``"blackwell"`` → Sv + split-beam angle detector + params : dict + Method-specific keyword arguments (see below). Unspecified keys use + the method's defaults. + + Method-specific arguments + ------------------------- + basic + var_name : str + Name of Sv variable (dB), e.g. ``"Sv"``. + channel : str + Channel identifier to select (e.g., ``"GPT 38 kHz ..."``). + threshold : float, default -50.0 + Sv threshold(s) in dB. If a single float is given, it is treated as the + lower bound and the upper bound is set to 10 dB above the lower + bound. If a 2-tuple `(tmin, tmax)` is provided, both the lower and + upper bounds are used directly. + offset_m : float, default 0.5 + Meters subtracted from the detected crossing. + bin_skip_from_surface : int, default 200 + Number of shallow range bins to ignore before searching. + + blackwell + var_name : str + Name of the Sv variable to use (e.g., ``"Sv"``). + channel : str + Channel identifier to select. + threshold : float | list | tuple, default -75 + Either a single Sv dB threshold (angle thresholds use defaults), + or a 3-tuple/list ``(tSv_dB, ttheta, tphi)`` after angle smoothing. + offset : float, default 0.3 + Meters subtracted from the detected bottom. + r0 : float, default 0 + Shallow bound (m) of the detection range. + r1 : float, default 500 + Deep bound (m) of the detection range. + wtheta : int, default 28 + Square smoothing window (pixels) for along-ship angle. + wphi : int, default 52 + Square smoothing window (pixels) for athwart-ship angle. + + Returns + ------- + xr.DataArray + 1-D bottom depth per ``ping_time`` (no ``channel`` dimension). The + output name and attributes are set by the called method. + + Raises + ------ + ValueError + If ``method`` is not supported. + + Examples + -------- + >>> detect_seafloor(ds, "basic", { + ... "var_name": "Sv", "channel": "GPT 38 kHz ...", + ... "threshold": -50, "offset_m": 0.5, "bin_skip_from_surface": 200 + ... }) + + >>> detect_seafloor(ds, "blackwell", { + ... "channel": "GPT 38 kHz ...", "threshold": (-75, 0.02, 0.02), + ... "offset": 0.3, "r0": 0, "r1": 600, "wtheta": 28, "wphi": 52 + ... }) + """ + if method not in METHODS_BOTTOM: + raise ValueError(f"Unsupported bottom detection method: {method}") + + return METHODS_BOTTOM[method](ds, **params) + + +# Registry of supported methods for shoal detection +METHODS_SHOAL = { + "echoview": shoal_echoview, + "weill": shoal_weill, +} + + +def detect_shoal( + ds: xr.Dataset, + method: str, + params: Dict, +) -> xr.DataArray: + """ + Detect shoals using the selected method and return a 2D boolean mask. + + Parameters + ---------- + ds : xr.Dataset + Sv dataset including ping_time and range_sample. + method : str + Name of the detection method to use (e.g., "echoview", "weill"). + params : dict + Parameters for the detection function (method-specific). + + Returns + ------- + xr.DataArray + 2D boolean DataArray of shoal mask (True = inside). + """ + if method not in METHODS_SHOAL: + raise ValueError(f"Unsupported shoal detection method: {method}") + + return METHODS_SHOAL[method](ds, **params) diff --git a/echopype/mask/seafloor_detection/bottom_basic.py b/echopype/mask/seafloor_detection/bottom_basic.py new file mode 100644 index 000000000..4872bc6ab --- /dev/null +++ b/echopype/mask/seafloor_detection/bottom_basic.py @@ -0,0 +1,100 @@ +import numpy as np +import xarray as xr + +from echopype.mask.seafloor_detection.utils import ( + _check_inputs, + _validate_threshold, +) + + +def bottom_basic( + ds: xr.Dataset, + var_name: str, + channel: str, + threshold: float = -50.0, + offset_m: float = 0.5, + bin_skip_from_surface: int = 200, +) -> xr.DataArray: + """ + Simple threshold-based seafloor detection returning a 1-D bottom line (depth). + + Summary + ------- + For the selected `channel`, the algorithm skips the top `bin_skip_from_surface` + range bins, then finds (per ping) the first range sample where `Sv` is within a + user-defined dB interval. That depth (looked up from `depth`) minus `offset_m` + is returned as the bottom. + + Parameters + ---------- + ds : xr.Dataset + Dataset containing: + • `var_name` (Sv in dB), typically with dims + (`channel`, `ping_time`, `range_sample`); + • a vertical coordinate (e.g., `depth`) aligned with `range_sample`. + var_name : str + Name of the Sv variable to use (e.g., `"Sv"`). + channel : str + Channel identifier to process (must match an entry in `ds['channel']`). + threshold : float or tuple(float, float), default -50.0 + Sv threshold(s) in dB. If a single float is given, it is treated as the + lower bound and the upper bound is set to 10 dB above the lower + bound. If a 2-tuple `(tmin, tmax)` is provided, both the lower + and upper bounds are used directly. + offset_m : float, default 0.5 + Meters subtracted from the detected crossing to place the bottom slightly + above the echo maximum. + bin_skip_from_surface : int, default 200 + Number of shallow range bins to ignore before searching (index units, + not meters). + + Returns + ------- + xr.DataArray + 1-D bottom depth per `ping_time` (no `channel` dimension) with attributes: + `detector='basic'`, `threshold_min`, `threshold_max`, `offset_m`, + `bin_skip_from_surface`, and `channel`. + + Notes + ----- + * Depth lookup uses `depth.isel(ping_time=0)` as the reference vector. + * If no sample meets the threshold in a ping, `argmax` on a False-only + mask returns 0 (i.e., the first searched bin after the skipped region). + """ + + Sv_sel, depth_sel = _check_inputs(ds, var_name, channel) + tmin, tmax = _validate_threshold(threshold) + + depth_ref = depth_sel.isel(ping_time=0) + Sv_sliced = Sv_sel.isel(range_sample=slice(bin_skip_from_surface, None)) + cond = (Sv_sliced > tmin) & (Sv_sliced < tmax) + + # Find index of first match along range_sample + idx = cond.argmax(dim="range_sample") + bin_skip_from_surface # add back skipped samples + + # Map index to depth using depth_ref + bottom_depth = xr.apply_ufunc( + lambda i: depth_ref[int(i)] if np.isfinite(i) else np.nan, + idx.astype(float), + vectorize=True, + dask="parallelized", + output_dtypes=[float], + ) - float( + offset_m + ) # we remove the offset + + # Return 1D DataArray with attributes + return xr.DataArray( + bottom_depth.data, + dims=["ping_time"], + coords={"ping_time": ds["ping_time"]}, + name="bottom_depth", + attrs={ + "detector": "basic", + "threshold_min": float(tmin), + "threshold_max": float(tmax), + "offset_m": float(offset_m), + "bin_skip_from_surface": int(bin_skip_from_surface), + "channel": str(channel), + }, + ) diff --git a/echopype/mask/seafloor_detection/bottom_blackwell.py b/echopype/mask/seafloor_detection/bottom_blackwell.py new file mode 100644 index 000000000..953bd78ef --- /dev/null +++ b/echopype/mask/seafloor_detection/bottom_blackwell.py @@ -0,0 +1,173 @@ +import numpy as np +import scipy.ndimage as ndima +import xarray as xr +from scipy.signal import convolve2d + +from echopype.mask.seafloor_detection.utils import _check_inputs, _parse_blackwell_thresholds +from echopype.utils.compute import _lin2log, _log2lin + + +def bottom_blackwell( + ds: xr.Dataset, + var_name: str, + channel: str, + threshold: float | list | tuple = -75, + offset: float = 0.3, + r0: float = 0, + r1: float = 500, + wtheta: int = 28, + wphi: int = 52, +) -> xr.DataArray: + """ + Shoal detector modified from the "blackwell" function + in `mask_seabed.py`, originally written by + Alejandro ARIZA for the Echopy library © 2020. + + Based on: "Blackwell et al (2019), Aliased seabed detection in fisheries acoustic + data". (https://arxiv.org/abs/1904.10736) + + Overview + --------- + 1) Range crop: restrict processing to r ∈ [r0, r1]. + + 2) Angle smoothing: convolve along-ship (θ) and athwart-ship (ϕ) angles + with square kernels of size wtheta and wphi (mean filters). + + 3) Angle activity mask: flag pixels where the smoothed angles’ squared values + exceed thresholds (θ > ttheta or ϕ > tphi), then take the union. + 4) Adaptive Sv threshold: compute the median Sv over the angle mask in linear + units (then convert to dB), set it to tSv if lower, and threshold Sv > Sv_median. + + 5) Connected components: label Sv above threshold patches and keep only those + intersecting the angle mask. + + 6) Bottom pick: for each ping, take the shallowest range of the kept patch and + subtract an `offset` (m) to place the bottom line slightly above it. + + Parameters + ---------- + ds : xr.Dataset + Dataset containing: + • ``var_name`` (Sv in dB) with dims typically + (``channel``, ``ping_time``, ``range_sample``), + • ``angle_alongship`` and ``angle_athwartship`` with compatible dims, + • a vertical coordinate (e.g., ``depth``) aligned with ``range_sample``. + var_name : str + Name of the Sv variable to use (e.g., ``"Sv"``). + channel : str + Channel identifier to process (must match an entry in ``ds['channel']``). + threshold : float | list | tuple, default -75 + Either a single Sv threshold in dB (angle thresholds use defaults), or a + 3-tuple/list ``(tSv_dB, ttheta, tphi)`` where ``ttheta`` and ``tphi`` are + post-smoothing angle activity thresholds (same units as the squared, + smoothed angles in your pipeline). + offset : float, default 0.3 + Meters subtracted from the detected range to place the bottom line slightly + above the echo maximum. + r0, r1 : float, default 0, 500 + Shallow and deep bounds (meters) of the search interval. + wtheta, wphi : int, default 28, 52 + Side length (pixels) of the square smoothing windows for the along-ship + and athwart-ship angle fields. + + Returns + ------- + xr.DataArray + 1-D bottom depth per ``ping_time`` with attributes: + ``detector='blackwell'``, ``threshold_Sv``, ``threshold_angle_major``, + ``threshold_angle_minor``, ``offset_m``, and ``channel``. + """ + + # Validate input variables and structure + Sv_sel, depth_sel = _check_inputs( + ds, + var_name=var_name, + channel=channel, + required_vars=["angle_alongship", "angle_athwartship"], + ) + + # parse thresholds + tSv, ttheta, tphi = _parse_blackwell_thresholds(threshold) + + # Direct selection without using ch_sel + theta = ds["angle_alongship"].sel(channel=channel) + phi = ds["angle_athwartship"].sel(channel=channel) + + # to match with blackwell echopy format + Sv_sel = Sv_sel.transpose("range_sample", "ping_time") + theta = theta.transpose("range_sample", "ping_time") + phi = phi.transpose("range_sample", "ping_time") + + ping_time = Sv_sel.coords["ping_time"] + r = depth_sel.isel(ping_time=0).values + + # Define core detection range + r0_idx = np.nanargmin(abs(r - r0)) + r1_idx = np.nanargmin(abs(r - r1)) + 1 + + Svchunk = Sv_sel[r0_idx:r1_idx, :] + thetachunk = theta[r0_idx:r1_idx, :] + phichunk = phi[r0_idx:r1_idx, :] + + # Build angle masks + ktheta = np.ones((wtheta, wtheta)) / wtheta**2 + kphi = np.ones((wphi, wphi)) / wphi**2 + + # Angle masks + thetamaskchunk = convolve2d(thetachunk, ktheta, "same", boundary="symm") ** 2 > ttheta + phimaskchunk = convolve2d(phichunk, kphi, "same", boundary="symm") ** 2 > tphi + anglemaskchunk = thetamaskchunk | phimaskchunk + + # Apply Blackwell algorithm + if anglemaskchunk.any(): + + Svmedian_anglemasked = float( + _lin2log(np.nanmedian(_log2lin(Svchunk.values[anglemaskchunk]))) + ) + + if np.isnan(Svmedian_anglemasked): + Svmedian_anglemasked = np.inf + if Svmedian_anglemasked < tSv: + Svmedian_anglemasked = tSv + + Svmaskchunk = Svchunk > Svmedian_anglemasked + + # Connected components + items = ndima.label(Svmaskchunk, ndima.generate_binary_structure(2, 2))[0] + intercepted = list(set(items[anglemaskchunk])) + if 0 in intercepted: + intercepted.remove(0) + + # Combine intercepted items + maskchunk = np.zeros(Svchunk.shape, dtype=bool) + for i in intercepted: + maskchunk |= items == i + + # Add padding + above = np.zeros((r0_idx, maskchunk.shape[1]), dtype=bool) + below = np.zeros((len(r) - r1_idx, maskchunk.shape[1]), dtype=bool) + + mask = np.r_[above, maskchunk, below] + + else: + mask = np.zeros_like(Sv_sel, dtype=bool) + + # Bottom detection from mask - offset + bottom_sample_idx = mask.argmax(axis=0) + bottom_depth = r[bottom_sample_idx] - offset + + # Return 1D DataArray with attributes + return xr.DataArray( + bottom_depth, + dims=["ping_time"], + coords={"ping_time": ping_time}, + name="bottom_depth", + attrs={ + "detector": "blackwell", + "threshold_Sv": float(tSv), + "threshold_angle_major": float(ttheta), + "threshold_angle_minor": float(tphi), + "offset_m": float(offset), + "channel": channel, + }, + ) diff --git a/echopype/mask/seafloor_detection/utils.py b/echopype/mask/seafloor_detection/utils.py new file mode 100644 index 000000000..e5e831f95 --- /dev/null +++ b/echopype/mask/seafloor_detection/utils.py @@ -0,0 +1,64 @@ +import xarray as xr + + +def _check_inputs(ds: xr.Dataset, var_name: str, channel: str, required_vars: list[str] = None): + """Validate dataset and select the reference channel for bottom detection.""" + if var_name not in ds: + raise KeyError(f"{var_name!r} not found in dataset") + if "depth" not in ds: + raise KeyError("'depth' variable not found in dataset") + if "channel" not in ds.coords: + raise ValueError("Dataset must have 'channel' coordinate") + + required_vars = required_vars or [] + for var in required_vars: + if var not in ds: + raise KeyError(f"Required variable {var!r} not found in dataset") + + Sv_all = ds[var_name] + depth_all = ds["depth"] + Sv_sel = Sv_all.sel(channel=channel) + depth_sel = depth_all.sel(channel=channel) + + # Ensure uniform depth grid + depth_ref = depth_sel.isel(ping_time=0) + is_uniform = (abs(depth_sel - depth_ref).max(dim="range_sample") < 1e-16).all() + if not bool(is_uniform): + raise ValueError("Depth grid varies across ping_time for the selected channel.") + + return Sv_sel, depth_sel + + +def _validate_threshold(threshold): + """Ensure threshold is a valid tuple (tmin, tmax).""" + if isinstance(threshold, (int, float)): + tmin, tmax = float(threshold), float(threshold) + 10.0 + else: + tmin, tmax = map(float, threshold) + if tmax <= tmin: + raise ValueError("threshold upper bound must be > lower bound") + return tmin, tmax + + +def _parse_blackwell_thresholds(threshold): + """ + Parse threshold for Blackwell detection. + + Returns + ------- + tuple : (tSv, ttheta, tphi) + Thresholds for Sv (dB), angle_major, angle_minor. + """ + if isinstance(threshold, (list, tuple)): + if len(threshold) == 3: + tSv, ttheta, tphi = threshold + elif len(threshold) == 2: + tSv, ttheta, tphi = threshold[0], 702, 282 + else: + raise ValueError("`threshold` must have 1, 2, or 3 values") + elif isinstance(threshold, (int, float)): + tSv, ttheta, tphi = threshold, 702, 282 + else: + raise TypeError("`threshold` must be float or tuple/list of 1–3 floats") + + return float(tSv), float(ttheta), float(tphi) diff --git a/echopype/mask/shoal_detection/shoal_echoview.py b/echopype/mask/shoal_detection/shoal_echoview.py new file mode 100644 index 000000000..65309cc1b --- /dev/null +++ b/echopype/mask/shoal_detection/shoal_echoview.py @@ -0,0 +1,144 @@ +import numpy as np +import pandas as pd +import scipy.ndimage as ndima +import xarray as xr + + +def shoal_echoview( + ds: xr.Dataset, + var_name: str, + channel: str | None, + idim: np.ndarray, + jdim: np.ndarray, + thr: float = -70.0, + mincan: tuple[float, float] = (3.0, 10.0), + maxlink: tuple[float, float] = (3.0, 15.0), + minsho: tuple[float, float] = (3.0, 15.0), +) -> xr.DataArray: + """ + Shoal detector modified from the "echoview" function + in `mask_shoals.py`, originally written by + Alejandro ARIZA for the Echopy library © 2020. + + Overview + --------- + Perform shoal detection on a Sv matrix using Echoview-like algorithm. + + 1) Threshold (candidates): mark samples where Sv > `thr`. + + 2) Minimum candidate size: remove connected components whose + height < `mincan[0]` (in idim units) or width < `mincan[1]` (in jdim units). + + 3) Linking: for each remaining component, search a surrounding box expanded + by `maxlink` (height, width) and link neighbouring components by assigning + them the same label. + + 4) Minimum shoal size: after linking, remove shoals whose height < `minsho[0]` + or width < `minsho[1]`. + + Parameters + ---------- + ds : xr.Dataset + Dataset containing `var_name` (Sv in dB). + var_name : str + Name of the Sv variable in `ds`. + channel : str | None + If a "channel" dimension exists, select this channel. + idim : np.ndarray + Vertical axis coordinates (length = n_rows + 1). + jdim : np.ndarray + Horizontal axis coordinates (length = n_cols + 1). + thr : float + Threshold in dB for initial detection (detect values > `thr`). + mincan : tuple[float, float] + Minimum candidate size (height, width). + maxlink : tuple[float, float] + Maximum linking distance (height, width). + minsho : tuple[float, float] + Minimum shoal size (height, width) after linking. + + Returns + ------- + xr.DataArray + Boolean mask of detected shoals (dims: "ping_time", "range_sample"; True = detected). + """ + + # Validate variable + if var_name not in ds: + raise ValueError(f"Variable '{var_name}' not found in dataset") + + var = ds[var_name] + + if "channel" in var.dims: + if channel is None: + raise ValueError("Please specify channel for multi-channel data") + var = var.sel(channel=channel) + + if np.isnan(idim).any() or np.isnan(jdim).any(): + raise ValueError("idim and jdim must not contain NaN") + + Sv = var.transpose("range_sample", "ping_time").values # (range, ping) + + # 1. Thresholding + mask = np.ma.masked_greater(Sv, thr).mask + if isinstance(mask, np.bool_): # scalar + mask = np.zeros_like(Sv, dtype=bool) + + # 2. Remove candidates below mincan size + candidateslabeled = ndima.label(mask, np.ones((3, 3)))[0] + candidateslabels = pd.factorize(candidateslabeled[candidateslabeled != 0])[1] + for cl in candidateslabels: + candidate = candidateslabeled == cl + idx = np.where(candidate)[0] + jdx = np.where(candidate)[1] + height = idim[max(idx + 1)] - idim[min(idx)] + width = jdim[max(jdx + 1)] - jdim[min(jdx)] + if (height < mincan[0]) or (width < mincan[1]): + mask[idx, jdx] = False + + # 3. Linking neighbours + linked = np.zeros(mask.shape, dtype=int) + shoalslabeled = ndima.label(mask, np.ones((3, 3)))[0] + shoalslabels = pd.factorize(shoalslabeled[shoalslabeled != 0])[1] + for fl in shoalslabels: + shoal = shoalslabeled == fl + i0, i1 = np.min(np.where(shoal)[0]), np.max(np.where(shoal)[0]) + j0, j1 = np.min(np.where(shoal)[1]), np.max(np.where(shoal)[1]) + i00 = np.nanargmin(abs(idim - (idim[i0] - (maxlink[0] + 1)))) + i11 = np.nanargmin(abs(idim - (idim[i1] + (maxlink[0] + 1)))) + 1 + j00 = np.nanargmin(abs(jdim - (jdim[j0] - (maxlink[1] + 1)))) + j11 = np.nanargmin(abs(jdim - (jdim[j1] + (maxlink[1] + 1)))) + 1 + + around = np.zeros_like(mask, dtype=bool) + around[i00:i11, j00:j11] = True + neighbours = around & mask + neighbourlabels = pd.factorize(shoalslabeled[neighbours])[1] + neighbourlabels = neighbourlabels[neighbourlabels != 0] + neighbours = np.isin(shoalslabeled, neighbourlabels) + + if (pd.factorize(linked[neighbours])[1] == 0).all(): + linked[neighbours] = np.max(linked) + 1 + else: + formerlabels = pd.factorize(linked[neighbours])[1] + minlabel = np.min(formerlabels[formerlabels != 0]) + linked[neighbours] = minlabel + for fl in formerlabels[formerlabels != 0]: + linked[linked == fl] = minlabel + + # 4. Remove linked shoals smaller than minsho + linkedlabels = pd.factorize(linked[linked != 0])[1] + for ll in linkedlabels: + shoal = linked == ll + idx, jdx = np.where(shoal) + height = idim[max(idx + 1)] - idim[min(idx)] + width = jdim[max(jdx + 1)] - jdim[min(jdx)] + if (height < minsho[0]) or (width < minsho[1]): + mask[idx, jdx] = False + + return xr.DataArray( + mask.T.astype(bool), + dims=["ping_time", "range_sample"], + coords={"ping_time": ds["ping_time"], "range_sample": ds["range_sample"]}, + name="shoal_mask", + attrs={"description": f"Shoal mask using Echoview algorithm on {var_name}"}, + ) diff --git a/echopype/mask/shoal_detection/shoal_weill.py b/echopype/mask/shoal_detection/shoal_weill.py new file mode 100644 index 000000000..5429000c2 --- /dev/null +++ b/echopype/mask/shoal_detection/shoal_weill.py @@ -0,0 +1,160 @@ +import numpy as np +import scipy.ndimage as ndi +import xarray as xr + + +def shoal_weill( + ds: xr.Dataset, + var_name: str, + channel: str | None = None, + thr: float = -70.0, + maxvgap: int = 5, + maxhgap: int = 0, + minvlen: int = 0, + minhlen: int = 0, +) -> xr.DataArray: + """ + Shoal detector modified from the "weill" function + in `mask_shoals.py`, originally written by + Alejandro ARIZA for the Echopy library © 2020. + + Detects and masks shoals following the algorithm described in: + + "Weill et al. (1993): MOVIES-B — an acoustic detection description + software. Application to shoal species' classification". + + Overview + --------- + Groups contiguous regions of Sv above a given threshold as a single shoal, + following the contiguity rules of Weill et al. (1993) in the following steps: + + 1) Threshold: mark samples where Sv > `thr`. + + 2) Vertical contiguity: within each ping, fill short gaps (≤ `maxvgap` range samples), + ignoring gaps that touch the top/bottom boundaries. + + 3) Horizontal contiguity: at each range, fill short gaps across pings (≤ `maxhgap` pings), + ignoring gaps that touch the left/right boundaries. + + 4) Size filter: remove features whose vertical length < `minvlen` (range samples) or + horizontal length < `minhlen` (pings). + + Parameters + ---------- + ds : xr.Dataset + Dataset containing `var_name` with Sv in dB. + var_name : str + Name of the Sv variable in `ds`. + channel : str | None + If a "channel" dimension exists, select this channel. + thr : float + Threshold in dB (detect values > `thr`). + maxvgap : int + Max vertical gap (in range samples) to fill. + maxhgap : int + Max horizontal gap (in pings) to fill. + minvlen : int + Minimum vertical length (in range samples) to keep a feature. + minhlen : int + Minimum horizontal length (in pings) to keep a feature. + + Returns + ------- + xr.DataArray + Boolean mask with dims ("ping_time", "range_sample") where True = detected. + """ + if var_name not in ds: + raise ValueError(f"Variable '{var_name}' not found in dataset") + + var = ds[var_name] + + # If multi-channel, select one + if "channel" in var.dims: + if channel is None: + raise ValueError("Please specify 'channel' for multi-channel data.") + var = var.sel(channel=channel) + + # Ensure we have the two core dims + if not {"ping_time", "range_sample"}.issubset(set(var.dims)): + raise ValueError( + f"'{var_name}' must have dims including 'ping_time' and 'range_sample', " + f"got {tuple(var.dims)}" + ) + + # Arrange as (range, ping) for processing, similar to echopy + Sv = var.transpose("range_sample", "ping_time").values + + # --- 1) Thresholding: keep (mask=True) where Sv > thr + mask = np.ma.masked_greater(Sv, thr).mask + if np.isscalar(mask): + mask = np.zeros_like(Sv, dtype=bool) + + n_range, n_ping = mask.shape + + # --- 2) Fill short vertical gaps per ping + # Work on each column (over range axis) + for jdx in range(n_ping): + col = mask[:, jdx] + # Label False regions (gaps) within the True mask + labelled = ndi.label(~col)[0] + # If all False or all True, nothing to do + if (labelled == 0).all() or (labelled == 1).all(): + continue + for lab in range(1, labelled.max() + 1): + gap = labelled == lab + gap_size = int(gap.sum()) + if gap_size <= maxvgap: + idx = np.where(gap)[0] + # Do not fill if the gap touches top/bottom boundary + if 0 in idx or (n_range - 1) in idx: + continue + mask[idx, jdx] = True + + # --- 3) Fill short horizontal gaps per depth + # Work on each row (over ping axis) + for idx in range(n_range): + row = mask[idx, :] + labelled = ndi.label(~row)[0] + if (labelled == 0).all() or (labelled == 1).all(): + continue + for lab in range(1, labelled.max() + 1): + gap = labelled == lab + gap_size = int(gap.sum()) + if gap_size <= maxhgap: + jdxs = np.where(gap)[0] + # Do not fill if the gap touches left/right boundary + if 0 in jdxs or (n_ping - 1) in jdxs: + continue + mask[idx, jdxs] = True + + # --- 4) Remove features smaller than (minvlen, minhlen) + # Label True regions and filter by size in (range, ping) coordinates + features = ndi.label(mask)[0] + if features.max() > 0: + for lab in range(1, features.max() + 1): + feat = features == lab + if not feat.any(): + continue + ii, jj = np.where(feat) + vlen = int(ii.max() - ii.min() + 1) # vertical length in samples + hlen = int(jj.max() - jj.min() + 1) # horizontal length in pings + if (vlen < minvlen) or (hlen < minhlen): + mask[ii, jj] = False + + # Return as (ping_time, range_sample) to match echopype convention + out = xr.DataArray( + mask.T.astype(bool), + dims=("ping_time", "range_sample"), + coords={"ping_time": ds["ping_time"], "range_sample": ds["range_sample"]}, + name="shoal_mask_weill", + attrs={ + "description": f"Weill-style threshold+gap-fill mask on '{var_name}'", + "threshold_dB": float(thr), + "maxvgap": int(maxvgap), + "maxhgap": int(maxhgap), + "minvlen": int(minvlen), + "minhlen": int(minhlen), + **({"channel": str(channel)} if channel is not None else {}), + }, + ) + return out diff --git a/echopype/test_data/README.md b/echopype/test_data/README.md deleted file mode 100644 index d3295604e..000000000 --- a/echopype/test_data/README.md +++ /dev/null @@ -1,31 +0,0 @@ -## Files Used in Tests -Most of these files are stored on Git LFS but the ones that aren't (due to file size) can be found on the echopype shared drive. - - -### EK80 - -- D20190822-T161221.raw: Contains channels that only record real power data -- D20170912-T234910.raw: Contains channels that only record complex power data -- Summer2018--D20180905-T033113.raw: Contains BB channels encoded in complex and CW channels encoded in power samples (reduced from 300 MB to 3.8 MB in test data updates). -- Summer2018 (3 files): Contains channels with complex power data as well as channels with real power data. They can be combined. -- 2019118 group2survey-D20191214-T081342.raw: Contains 6 channels but only 2 of those channels collect ping data -- D20200528-T125932.raw: Data collected from WBT mini (instead of WBT), from @emlynjdavies -- Green2.Survey2.FM.short.slow.-D20191004-T211557.raw: Contains 2-in-1 transducer, from @FletcherFT (reduced from 104.9 MB to 765 KB in test data updates) -- raw4-D20220514-T172704.raw: Contains RAW4 datagram, 1 channel only, from @cornejotux -- D20210330-T123857.raw: do not contain filter coefficients - - -### EA640 -- 0001a-D20200321-T032026.raw: Data of identical format to standard EK80 files, but with a different NME1 datagram (instead of NME0). - - -### EK60 -- DY1801_EK60-D20180211-T164025.raw: Standard test with constant ranges across ping times -- Winter2017-D20170115-T150122.raw: Contains a change of recording length in the middle of the file -- 2015843-D20151023-T190636.raw: Not used in tests but contains ranges are not constant across ping times -- SH1701_consecutive_files_w_range_change: Not used in tests. [Folder](https://drive.google.com/drive/u/1/folders/1PaDtL-xnG5EK3N3P1kGlXa5ub16Yic0f) on shared drive that contains sequential files with ranges that are not constant across ping times. -- NBP_B050N-D20180118-T090228.raw: split-beam setup without angle data - - -### AZFP -- 17082117.01A: Standard test. Used with 17041823.XML. diff --git a/echopype/testing.py b/echopype/testing.py index 75a19d4e7..44259ed5b 100644 --- a/echopype/testing.py +++ b/echopype/testing.py @@ -1,5 +1,4 @@ from collections import defaultdict -from pathlib import Path import numpy as np import pandas as pd @@ -7,9 +6,6 @@ from .utils.compute import _lin2log, _log2lin -HERE = Path(__file__).parent.absolute() -TEST_DATA_FOLDER = HERE / "test_data" - # Data length for each data type _DATA_LEN = { "power": 1, @@ -24,7 +20,7 @@ def _gen_ping_time(ping_time_len, ping_time_interval, ping_time_jitter_max_ms=0) jitter = ( np.random.randint(ping_time_jitter_max_ms, size=ping_time_len) / 1000 ) # convert to seconds - ping_time = pd.to_datetime(ping_time.astype(int) / 1e9 + jitter, unit="s") + ping_time = pd.to_datetime(ping_time.astype("int64") / 1e9 + jitter, unit="s") return ping_time diff --git a/echopype/tests/calibrate/test_calibrate.py b/echopype/tests/calibrate/test_calibrate.py index 3e5f28c15..98aa333b3 100644 --- a/echopype/tests/calibrate/test_calibrate.py +++ b/echopype/tests/calibrate/test_calibrate.py @@ -18,6 +18,9 @@ def azfp_path(test_path): def ek60_path(test_path): return test_path['EK60'] +@pytest.fixture +def ek60_cal_chunks_path(test_path): + return test_path['EK60_CAL_CHUNKS'] @pytest.fixture def ek80_path(test_path): @@ -270,7 +273,7 @@ def test_compute_Sv_ek80_CW_complex_BB_complex(ek80_cal_path, ek80_path): @pytest.mark.integration -def test_compute_Sv_combined_ed_ping_time_extend_past_time1(): +def test_compute_Sv_combined_ed_ping_time_extend_past_time1(ek80_path): """ Test computing combined Echodata object when ping time dimension in Beam group extends past time1 dimension in Environment group. @@ -280,8 +283,8 @@ def test_compute_Sv_combined_ed_ping_time_extend_past_time1(): """ # Parse RAW files and combine Echodata objects raw_list = [ - "echopype/test_data/ek80/pifsc_saildrone/SD_TPOS2023_v03-Phase0-D20230530-T001150-0.raw", - "echopype/test_data/ek80/pifsc_saildrone/SD_TPOS2023_v03-Phase0-D20230530-T002350-0.raw", + ek80_path / "pifsc_saildrone/SD_TPOS2023_v03-Phase0-D20230530-T001150-0.raw", + ek80_path / "pifsc_saildrone/SD_TPOS2023_v03-Phase0-D20230530-T002350-0.raw", ] ed_list = [] for raw_file in raw_list: @@ -321,11 +324,11 @@ def test_compute_Sv_combined_ed_ping_time_extend_past_time1(): @pytest.mark.parametrize( "raw_path, sonar_model, xml_path, waveform_mode, encode_mode", [ - ("azfp/17031001.01A", "AZFP", "azfp/17030815.XML", None, None), - ("ek60/DY1801_EK60-D20180211-T164025.raw", "EK60", None, None, None), - ("ek80/D20170912-T234910.raw", "EK80", None, "BB", "complex"), - ("ek80/D20230804-T083032.raw", "EK80", None, "CW", "complex"), - ("ek80/Summer2018--D20180905-T033113.raw", "EK80", None, "CW", "power") + ("17031001.01A", "AZFP", "17030815.XML", None, None), + ("DY1801_EK60-D20180211-T164025.raw", "EK60", None, None, None), + ("D20170912-T234910.raw", "EK80", None, "BB", "complex"), + ("D20230804-T083032.raw", "EK80", None, "CW", "complex"), + ("Summer2018--D20180905-T033113.raw", "EK80", None, "CW", "power") ] ) def test_check_echodata_backscatter_size( @@ -334,15 +337,29 @@ def test_check_echodata_backscatter_size( xml_path, waveform_mode, encode_mode, - caplog + caplog, + azfp_path, + ek60_path, + ek80_path ): """Tests for _check_echodata_backscatter_size warning.""" # Parse Echodata Object - ed = ep.open_raw( - raw_file=f"echopype/test_data/{raw_path}", - sonar_model=sonar_model, - xml_path=f"echopype/test_data/{xml_path}", - ) + if sonar_model == "AZFP": + ed = ep.open_raw( + raw_file=azfp_path / raw_path, + sonar_model=sonar_model, + xml_path=azfp_path / xml_path, + ) + elif sonar_model == "EK60": + ed = ep.open_raw( + raw_file=ek60_path / raw_path, + sonar_model=sonar_model, + ) + else: # EK80 + ed = ep.open_raw( + raw_file=ek80_path / raw_path, + sonar_model=sonar_model, + ) # Compute environment parameters if AZFP env_params = None @@ -420,10 +437,10 @@ def test_check_echodata_backscatter_size( @pytest.mark.integration -def test_fm_equals_bb(): +def test_fm_equals_bb(ek80_path): """Check that waveform_mode='BB' and waveform_mode='FM' result in the same Sv/TS.""" # Open Raw and Compute both Sv and both TS - ed = ep.open_raw("echopype/test_data/ek80/D20170912-T234910.raw", sonar_model = "EK80") + ed = ep.open_raw(ek80_path / "D20170912-T234910.raw", sonar_model = "EK80") ds_Sv_bb = ep.calibrate.compute_Sv(ed, waveform_mode="BB", encode_mode="complex") ds_Sv_fm = ep.calibrate.compute_Sv(ed, waveform_mode="FM", encode_mode="complex") ds_TS_bb = ep.calibrate.compute_TS(ed, waveform_mode="BB", encode_mode="complex") @@ -435,16 +452,15 @@ def test_fm_equals_bb(): @pytest.mark.integration -def test_calibrate_ek60_chunks(): +def test_calibrate_ek60_chunks(ek60_cal_chunks_path): """ Tests that the correct chunks are generated in the intermediate 'get_vend_cal_params_power' output and in the final 'compute_Sv' output. """ # Lazy-load the Zarr stores - test_dir = "echopype/test_data/ek60_calibrate_chunks" - ed = ep.open_converted(f"{test_dir}/ed_ek60_chunk_test.zarr", chunks={}) - beam = xr.open_zarr(f"{test_dir}/beam_ek60_chunk_test.zarr", chunks={}) - vend = xr.open_zarr(f"{test_dir}/vend_ek60_chunk_test.zarr", chunks={}) + ed = ep.open_converted(ek60_cal_chunks_path / "ed_ek60_chunk_test.zarr", chunks={}) + beam = xr.open_zarr(ek60_cal_chunks_path / "beam_ek60_chunk_test.zarr", chunks={}) + vend = xr.open_zarr(ek60_cal_chunks_path / "vend_ek60_chunk_test.zarr", chunks={}) # Get gain correction parameters stored in the Vendor_specific group da_param = get_vend_cal_params_power(beam=beam, vend=vend, param="gain_correction") diff --git a/echopype/tests/calibrate/test_calibrate_ek80.py b/echopype/tests/calibrate/test_calibrate_ek80.py index 75423bc06..484b1cc9e 100644 --- a/echopype/tests/calibrate/test_calibrate_ek80.py +++ b/echopype/tests/calibrate/test_calibrate_ek80.py @@ -21,6 +21,9 @@ def ek80_cal_path(test_path): def ek80_ext_path(test_path): return test_path["EK80_EXT"] +@pytest.fixture +def ek80_multiplex_path(test_path): + return test_path["EK80_MULTIPLEX"] def test_ek80_transmit_chirp(ek80_cal_path, ek80_ext_path): """ @@ -370,7 +373,7 @@ def test_ek80_BB_power_echoview(ek80_path): ev_vals = df_real.values[:, :] ep_vals = pc_mean.values.real[:, :] assert np.allclose(ev_vals[:, 69:], ep_vals[:, 69:], atol=1e-4) - assert np.allclose(ev_vals[:, 90:], ep_vals[:, 90:], atol=1e-5) + assert np.allclose(ev_vals[:, 90:], ep_vals[:, 90:], atol=1e-4) def test_ek80_CW_complex_Sv_receiver_sampling_freq(ek80_path): ek80_raw_path = str(ek80_path.joinpath("D20230804-T083032.raw")) @@ -417,9 +420,9 @@ def test_ek80_CW_complex_Sv_receiver_sampling_freq(ek80_path): ], ) @pytest.mark.integration -def test_ek80_BB_complex_multiplex_NaNs_and_non_NaNs(raw_data_path, target_channel_ping_pattern): +def test_ek80_BB_complex_multiplex_NaNs_and_non_NaNs(raw_data_path, target_channel_ping_pattern, ek80_multiplex_path): # Extract bb complex multiplex EK80 data - ed = ep.open_raw(f"echopype/test_data/ek80_bb_complex_multiplex/{raw_data_path}", sonar_model="EK80") + ed = ep.open_raw(ek80_multiplex_path / raw_data_path, sonar_model="EK80") # Compute Sv ds_Sv = ep.calibrate.compute_Sv(ed,waveform_mode='BB',encode_mode='complex') diff --git a/echopype/tests/calibrate/test_ecs_integration.py b/echopype/tests/calibrate/test_ecs_integration.py index 82e2d4b4f..46a6feab5 100644 --- a/echopype/tests/calibrate/test_ecs_integration.py +++ b/echopype/tests/calibrate/test_ecs_integration.py @@ -106,14 +106,15 @@ def test_ecs_intake_ek80_CW_power(ek80_path, ecs_path): # Check the final stored params (which are those used in calibration operations) # For those pulled from ECS for p_name in ["sound_speed", "temperature", "salinity", "pressure", "pH"]: - assert ds_Sv[p_name].identical(ecs_env_params[p_name]) + assert ds_Sv[p_name].equals(ecs_env_params[p_name]) + # replaced identical() with equals() because newer xarray versions check attrs for p_name in [ "sa_correction", "gain_correction", "equivalent_beam_angle", "beamwidth_alongship", "beamwidth_athwartship", "angle_offset_alongship", "angle_offset_athwartship", "angle_sensitivity_alongship", "angle_sensitivity_athwartship" ]: - assert ds_Sv[p_name].identical(ecs_cal_params[p_name]) + assert ds_Sv[p_name].equals(ecs_cal_params[p_name]) # For those computed from values in ECS file assert np.all(ds_Sv["sound_absorption"].values == assimilated_env_params["sound_absorption"].values) diff --git a/echopype/tests/clean/test_noise.py b/echopype/tests/clean/test_noise.py index e792df4fd..0a3620a29 100644 --- a/echopype/tests/clean/test_noise.py +++ b/echopype/tests/clean/test_noise.py @@ -13,6 +13,9 @@ ) from echopype.utils.compute import _lin2log, _log2lin +@pytest.fixture +def ek60_path(test_path): + return test_path["EK60"] @pytest.mark.unit def test_extract_dB(): @@ -43,11 +46,11 @@ def test_extract_dB(): ("echo_range"), ], ) -def test_mask_functions_with_no_vertical_range_variables(range_var): +def test_mask_functions_with_no_vertical_range_variables(range_var, ek60_path): """Test mask functions when no vertical range variables are in `ds_Sv`.""" # Open raw and calibrate ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw", + ek60_path / "from_echopy/JR230-D20091215-T121917.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -67,11 +70,11 @@ def test_mask_functions_with_no_vertical_range_variables(range_var): @pytest.mark.integration -def test_mask_functions_dimensions(): +def test_mask_functions_dimensions(ek60_path): """Test mask functions' output dimensions.""" # Open raw, calibrate, and add depth ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw", + ek60_path / "from_echopy/JR230-D20091215-T121917.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -91,11 +94,11 @@ def test_mask_functions_dimensions(): @pytest.mark.integration -def test_transient_mask_noise_func_error_and_warnings(caplog): +def test_transient_mask_noise_func_error_and_warnings(caplog, ek60_path): """Check if appropriate warnings and errors are raised for transient noise mask func input.""" # Open raw, calibrate, and add depth ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw", + ek60_path / "from_echopy/JR230-D20091215-T121917.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -155,14 +158,14 @@ def test_transient_mask_noise_func_error_and_warnings(caplog): (True, np.nanmedian), ], ) -def test_pool_Sv_values(chunk, func): +def test_pool_Sv_values(chunk, func, ek60_path): """ Manually check if the pooled Sv for transient noise masking contains the correct nan boundary and the correct bin aggregate values. """ # Open raw, calibrate, and add depth ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw", + ek60_path / "from_echopy/JR230-D20091215-T121917.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -258,11 +261,11 @@ def test_pool_Sv_values(chunk, func): (True, "nanmedian"), ], ) -def test_transient_noise_mask_values(chunk, func): +def test_transient_noise_mask_values(chunk, func, ek60_path): """Manually check if impulse noise mask removes transient noise values.""" # Open raw, calibrate, and add depth ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw", + ek60_path / "from_echopy/JR230-D20091215-T121917.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -336,7 +339,7 @@ def test_transient_noise_mask_values(chunk, func): (True, np.nanmedian), ], ) -def test_index_binning_pool_Sv_values(chunk, func): +def test_index_binning_pool_Sv_values(chunk, func, ek60_path): """ Manually check if the index binning pooled Sv for transient noise masking does the correct reflection computation. This is tested using `np.pad` to extend the Sv boundary @@ -344,7 +347,7 @@ def test_index_binning_pool_Sv_values(chunk, func): based off of this extended Sv. """ # Open raw, calibrate, and add depth - ed = ep.open_raw("echopype/test_data/ek60/from_echopy/JR161-D20061118-T010645.raw", sonar_model="EK60") + ed = ep.open_raw(ek60_path / "from_echopy/JR161-D20061118-T010645.raw", sonar_model="EK60") ds_Sv = ep.calibrate.compute_Sv(ed) ds_Sv = ep.consolidate.add_depth(ds_Sv) @@ -449,10 +452,10 @@ def test_index_binning_pool_Sv_values(chunk, func): (True, "nanmedian"), ], ) -def test_index_binning_transient_noise_mask_values(chunk, func): +def test_index_binning_transient_noise_mask_values(chunk, func, ek60_path): """Manually check if impulse noise mask removes transient noise values when using index binning.""" # Open raw, calibrate, and add depth - ed = ep.open_raw("echopype/test_data/ek60/from_echopy/JR161-D20061118-T010645.raw", sonar_model="EK60") + ed = ep.open_raw(ek60_path / "from_echopy/JR161-D20061118-T010645.raw", sonar_model="EK60") ds_Sv = ep.calibrate.compute_Sv(ed) ds_Sv = ep.consolidate.add_depth(ds_Sv) @@ -544,11 +547,11 @@ def test_index_binning_transient_noise_mask_values(chunk, func): (True), ], ) -def test_downsample_upsample_along_depth(chunk): +def test_downsample_upsample_along_depth(chunk, ek60_path): """Test downsample bins and upsample repeating values""" # Open raw, calibrate, and add depth ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw", + ek60_path / "from_echopy/JR230-D20091215-T121917.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -610,11 +613,11 @@ def test_downsample_upsample_along_depth(chunk): (True), ], ) -def test_index_binning_downsample_upsample_along_depth(chunk): +def test_index_binning_downsample_upsample_along_depth(chunk, ek60_path): """Test index binning downsampled-upsampled values.""" # Open raw, calibrate, and add depth ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw", + ek60_path / "from_echopy/JR230-D20091215-T121917.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -690,11 +693,11 @@ def test_index_binning_downsample_upsample_along_depth(chunk): (True, True), ], ) -def test_impulse_noise_mask_values(chunk, use_index_binning): +def test_impulse_noise_mask_values(chunk, use_index_binning, ek60_path): """Manually check if impulse noise mask removes impulse noise values.""" # Open raw, calibrate, and add depth ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw", + ek60_path / "from_echopy/JR230-D20091215-T121917.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -766,11 +769,11 @@ def test_impulse_noise_mask_values(chunk, use_index_binning): @pytest.mark.integration -def test_mask_attenuated_signal_limit_error(): +def test_mask_attenuated_signal_limit_error(ek60_path): """Test `mask_attenuated_signal` limit error.""" # Parse, calibrate, and add depth ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR161-D20061118-T010645.raw", + ek60_path / "from_echopy/JR161-D20061118-T010645.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -786,11 +789,11 @@ def test_mask_attenuated_signal_limit_error(): @pytest.mark.integration -def test_mask_attenuated_signal_outside_searching_range(): +def test_mask_attenuated_signal_outside_searching_range(ek60_path): """Test `mask_attenuated_signal` values errors.""" # Parse, calibrate, and add_depth ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR161-D20061118-T010645.raw", + ek60_path / "from_echopy/JR161-D20061118-T010645.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -820,11 +823,11 @@ def test_mask_attenuated_signal_outside_searching_range(): (True), ], ) -def test_mask_attenuated_signal_against_echopy(chunk): +def test_mask_attenuated_signal_against_echopy(chunk, ek60_path): """Test `attenuated_signal` to see if Echopype output matches echopy output mask.""" # Parse, calibrate, and add depth ed = ep.open_raw( - "echopype/test_data/ek60/from_echopy/JR161-D20061118-T010645.raw", + ek60_path / "from_echopy/JR161-D20061118-T010645.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) @@ -848,7 +851,7 @@ def test_mask_attenuated_signal_against_echopy(chunk): # Grab echopy attenuated signal mask echopy_attenuated_mask = xr.open_dataset( - "echopype/test_data/ek60/from_echopy/JR161-D20061118-T010645_echopy_attenuated_masks.zarr", + ek60_path / "from_echopy/JR161-D20061118-T010645_echopy_attenuated_masks.zarr", engine="zarr" ) @@ -860,11 +863,11 @@ def test_mask_attenuated_signal_against_echopy(chunk): @pytest.mark.unit -def test_estimate_background_noise_upsampling(): +def test_estimate_background_noise_upsampling(ek60_path): """Test for the correct upsampling behavior in `estimate_background_noise`""" # Open Raw and Calibrate ed = ep.open_raw( - "echopype/test_data/ek60/ncei-wcsd/Summer2017-D20170615-T190214.raw", + ek60_path / "ncei-wcsd/Summer2017-D20170615-T190214.raw", sonar_model="EK60" ) ds_Sv = ep.calibrate.compute_Sv(ed) diff --git a/echopype/tests/clean/test_transient_noise.py b/echopype/tests/clean/test_transient_noise.py new file mode 100644 index 000000000..6393c88c7 --- /dev/null +++ b/echopype/tests/clean/test_transient_noise.py @@ -0,0 +1,176 @@ +import numpy as np +import xarray as xr +import echopype as ep +import pytest + +# ---------- Fixtures + +@pytest.fixture(scope="module") +def ek60_path(test_path): + return test_path["EK60"] + +@pytest.fixture(scope="module") +def ds_small(ek60_path): + """Open raw, calibrate to Sv, add depth, and take a small deterministic slice.""" + ed = ep.open_raw( + ek60_path / "DY1801_EK60-D20180211-T164025.raw", + sonar_model="EK60", + ) + ds_Sv = ep.calibrate.compute_Sv(ed) + ds_Sv = ep.consolidate.add_depth(ds_Sv) + return ds_Sv + +# don t know if useful at the moment with code implementation +@pytest.fixture(params=[False, True], ids=["unchunked", "chunked"]) +def ds_small_chunked(ds_small, request): + """Parametrize chunking.""" + return ds_small.chunk("auto") if request.param else ds_small + + +# ---------- Dispatcher tests + +@pytest.mark.unit +def test_dispatcher_rejects_unsupported_method(ds_small): + with pytest.raises(ValueError, match="Unsupported transient noise removal method"): + ep.clean.detect_transient( + ds_small, method="not_a_method", params={} + ) + + +@pytest.mark.unit +@pytest.mark.parametrize("method,expected_name", [ + ("fielding", "fielding_mask_valid"), + ("matecho", "matecho_mask_valid"), +]) +def test_dispatcher_returns_named_boolean_mask(ds_small, method, expected_name): + params = { + # generic minimal args that both methods accept + "range_var": "depth", + } + mask = ep.clean.detect_transient(ds_small, method=method, params=params) + assert isinstance(mask, xr.DataArray) + assert mask.dtype == bool + assert mask.name == expected_name + # Dims and coords should match Sv exactly + Sv = ds_small["Sv"] + assert tuple(mask.dims) == tuple(Sv.dims) + for dim in Sv.dims: + assert Sv[dim].equals(mask[dim]) + + +# ---------- Fielding method tests + +@pytest.mark.integration +def test_fielding_dimensions_and_determinism(ds_small_chunked): + params = dict( + range_var="depth", + r0=900, r1=1000, n=30, thr=(3, 1), roff=20, jumps=5, maxts=-35, start=0, + ) + m1 = ep.clean.detect_transient(ds_small_chunked, "fielding", params) + m2 = ep.clean.detect_transient(ds_small_chunked, "fielding", params) + # Dims equal Sv + Sv = ds_small_chunked["Sv"] + assert tuple(m1.dims) == tuple(Sv.dims) + # Deterministic + xr.testing.assert_identical(m1, m2) + + +@pytest.mark.integration +def test_fielding_invalid_inputs_raise(ds_small): + # missing var_name + bad_ds = ds_small.drop_vars("Sv") + with pytest.raises(ValueError): + ep.clean.detect_transient(bad_ds, "fielding", dict(range_var="depth")) + # missing range_var + bad_ds2 = ds_small.drop_vars("depth") + with pytest.raises(ValueError): + ep.clean.detect_transient(bad_ds2, "fielding", dict(range_var="depth")) + + +# ---------- Matecho method tests + +@pytest.mark.integration +def test_matecho_dimensions_and_determinism(ds_small_chunked): + params = dict( + range_var="depth", + start_depth=220, + window_meter=450, + window_ping=100, + percentile=25, + delta_db=12, + extend_ping=0, + min_window=20, + ) + m1 = ep.clean.detect_transient(ds_small_chunked, "matecho", params) + m2 = ep.clean.detect_transient(ds_small_chunked, "matecho", params) + Sv = ds_small_chunked["Sv"] + assert tuple(m1.dims) == tuple(Sv.dims) + xr.testing.assert_identical(m1, m2) + + +@pytest.mark.integration +def test_matecho_threshold_monotonicity(ds_small): + """ + Increasing delta_db should make the detector more permissive: + i.e., fewer columns flagged as transient → more True (valid) in the mask. + """ + base_params = dict( + range_var="depth", + start_depth=220, + window_meter=450, + window_ping=100, + percentile=25, + extend_ping=0, + min_window=20, + ) + m_low = ep.clean.detect_transient(ds_small, "matecho", dict(delta_db=8, **base_params)) + m_high = ep.clean.detect_transient(ds_small, "matecho", dict(delta_db=16, **base_params)) + + # Masks are True=VALID. A higher threshold should yield >= number of True. + valid_low = np.count_nonzero(m_low.values) + valid_high = np.count_nonzero(m_high.values) + assert valid_high >= valid_low + + +@pytest.mark.integration +def test_matecho_bottom_var_optional(ds_small): + """ + Matecho accepts missing/NaN bottom and should still run. + If a (shallow) bottom is provided, it must not crash and must keep dims. + """ + # Run without bottom (default path in your wrapper) + m_nob = ep.clean.detect_transient( + ds_small, "matecho", + dict(range_var="depth", start_depth=220, window_meter=450, window_ping=50, delta_db=12) + ) + Sv = ds_small["Sv"] + assert tuple(m_nob.dims) == tuple(Sv.dims) + + # Provide a synthetic shallow bottom to exercise the code path + shallow = xr.DataArray( + np.full(ds_small.dims["ping_time"], 300.0), + dims=["ping_time"], + coords=[ds_small["ping_time"]], + name="bottom", + ) + ds_btm = ds_small.assign(bottom_var=shallow) + + # Your current wrapper ignores bottom_var argument; here we at least ensure it doesn't explode + m_btm = ep.clean.detect_transient( + ds_btm, "matecho", + dict(range_var="depth", start_depth=220, window_meter=450, window_ping=50, delta_db=12) + ) + assert tuple(m_btm.dims) == tuple(Sv.dims) + + +# ---------- Cross-method consistency + +@pytest.mark.integration +def test_methods_return_boolean_and_same_shape(ds_small_chunked): + params_fielding = dict(range_var="depth") + params_matecho = dict(range_var="depth") + mf = ep.clean.detect_transient(ds_small_chunked, "fielding", params_fielding) + mm = ep.clean.detect_transient(ds_small_chunked, "matecho", params_matecho) + Sv = ds_small_chunked["Sv"] + assert mf.dtype == bool and mm.dtype == bool + assert tuple(mf.dims) == tuple(Sv.dims) == tuple(mm.dims) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 8efd2d0d8..20dc6c36c 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -256,6 +256,22 @@ def test_compute_MVBS_invalid_range_var(ds_Sv_echo_range_regular, range_var): else: pass +@pytest.mark.integration +def test_compute_MVBS_w_dim_swapped(request): + """ + Test swapping dim from "channel" to "frequency_nominal" then computing MVBS + """ + ds_Sv = request.getfixturevalue("mock_Sv_dataset_regular") + ping_time_bin = "20s" + sv_swapped = ep.consolidate.swap_dims_channel_frequency(ds_Sv) + ds_MVBS = ep.commongrid.compute_MVBS(sv_swapped, ping_time_bin=ping_time_bin) + assert ds_MVBS is not None + + # Test to see if ping_time was resampled correctly + expected_ping_time = ( + ds_Sv["ping_time"].resample(ping_time=ping_time_bin, skipna=True).asfreq().indexes["ping_time"] + ) + assert np.array_equal(ds_MVBS.ping_time.data, expected_ping_time.values) @pytest.mark.integration def test_compute_MVBS(test_data_samples): @@ -316,10 +332,11 @@ def test_compute_MVBS_range_output(request, er_type): ds_MVBS = ep.commongrid.compute_MVBS(ds_Sv, range_bin="5m", ping_time_bin="10s") if er_type == "regular": + dt_ns = np.diff(ds_Sv["ping_time"][[0, -1]].values)[0] expected_len = ( ds_Sv["channel"].size, # channel - np.ceil(np.diff(ds_Sv["ping_time"][[0, -1]].astype(int)) / 1e9 / 10), # ping_time - np.ceil(ds_Sv["echo_range"].max() / 5), # depth + int(np.ceil(int(dt_ns) / 1e9 / 10)), # ping_time + int(np.ceil(ds_Sv["echo_range"].max() / 5)), # depth ) assert ds_MVBS["Sv"].shape == expected_len else: diff --git a/echopype/tests/conftest.py b/echopype/tests/conftest.py index 5e51c92cc..e895e4327 100644 --- a/echopype/tests/conftest.py +++ b/echopype/tests/conftest.py @@ -1,9 +1,139 @@ -"""``pytest`` configuration.""" +"""pytest configuration with minimal Pooch fallback for CI""" +import os import pytest +from pathlib import Path +if os.getenv("USE_POOCH") == "True": + import pooch -from echopype.testing import TEST_DATA_FOLDER + # Lock to the known-good assets release (can be overridden via env if needed) + ver = os.getenv("ECHOPYPE_DATA_VERSION", "v0.11.1a2") + base = os.getenv( + "ECHOPYPE_DATA_BASEURL", + "https://github.com/OSOceanAcoustics/echopype/releases/download/{version}/", + ) + cache_dir = pooch.os_cache("echopype") + + print( + "\n[echopype-ci] POOCH CONFIG\n" + f" USE_POOCH = {os.getenv('USE_POOCH')}\n" + f" ECHOPYPE_DATA_VERSION = {ver}\n" + f" ECHOPYPE_DATA_BASEURL = {base}\n" + f" pooch cache_dir = {cache_dir}\n", + flush=True, + ) + + bundles = [ + "ad2cp.zip", "azfp.zip", "azfp6.zip", "ea640.zip", "ecs.zip", "ek60.zip", + "ek60_calibrate_chunks.zip", "ek60_missing_channel_power.zip", "ek80.zip", + "ek80_bb_complex_multiplex.zip", "ek80_bb_with_calibration.zip", + "ek80_duplicate_ping_times.zip", "ek80_ext.zip", "ek80_invalid_env_datagrams.zip", + "ek80_missing_sound_velocity_profile.zip", "ek80_new.zip", "ek80_sequence.zip", + "es60.zip", "es70.zip", "es80.zip", "legacy_datatree.zip", + ] + + # v0.11.1a2 checksums (GitHub release assets) + registry = { + "ad2cp.zip": "sha256:78c634c7345991177b267c4cbb31f391990d2629b7f4a546da20d5126978b98a", + "azfp.zip": "sha256:fc75b48c81f266ce70d9db79a986fe8de4399c93bef35119acdc17a2d84aed49", + "azfp6.zip": "sha256:98228329333064fb4b44d3044296c79d58ac22f6d81f7f22cf770bacf0e882fd", + "ea640.zip": "sha256:49f70bd6f2355cb3c4c7a5b31fc00f7ae8c8a9ae888f0df1efe759032f9580df", + "ecs.zip": "sha256:dcc312baa1e9da4488f33bef625b1f86c8a92e3262e34fc90ccd0a4f90d1e313", + "ek60.zip": "sha256:66735de0ac584ec8a150b54b1a54024a92195f64036134ffdc9d472d7e155bb2", + "ek60_calibrate_chunks.zip": "sha256:bf435b1f7fc055f51afd55c4548713ba8e1eb0e919a0d74f4b9dd5f60b7fe327", + "ek60_missing_channel_power.zip": "sha256:f3851534cdc6ad3ae1d7c52a11cb279305d316d0086017a305be997d4011e20e", + "ek80.zip": "sha256:a114a8272e4c0e08c45c75241c50e3fd9e954f85791bb5eda25be98f6f782397", + "ek80_bb_complex_multiplex.zip": "sha256:9feaa90ce60036db8110ff14e732910786c6deff96836f2d6504ec1e9de57533", + "ek80_bb_with_calibration.zip": "sha256:53f018b6dae051cc86180e13cb3f28848750014dfcf84d97cf2191be2b164ccb", + "ek80_duplicate_ping_times.zip": "sha256:11a2dcb5cf113fa1bb03a6724524ac17bdb0db66cb018b0a3ca7cad87067f4bb", + "ek80_ext.zip": "sha256:79dd12b2d9e0399f88c98ab53490f5d0a8d8aed745650208000fcd947dbdd0de", + "ek80_invalid_env_datagrams.zip": "sha256:dece27d90f30d1a13b56d99350c3254e81622af3199fda0112d3b9e1d7db270c", + "ek80_missing_sound_velocity_profile.zip": "sha256:1635585026ae5c4ffdff09ca4d63aeff0b33471c5ee0e1b8a520f87469535852", + "ek80_new.zip": "sha256:f799cde453762c46ad03fee178c76cd9fbb00eec92a5d1038c32f6a9479b2e57", + "ek80_sequence.zip": "sha256:9d8fac39dd31f587d55b9978ba4d2b52bbc85daa85d320ef2ac34b3ae947bb1f", + "es60.zip": "sha256:a6c2a15c664ef8b6ac17cb36a28162c271fca361509cf43313038f1bdc9b6c7c", + "es70.zip": "sha256:a6b4f27f33f09bace26264de6984fdb4111a3a0337bc350c3c1d25c8b3effc7c", + "es80.zip": "sha256:b37ee01462f46efe055702c20be67d2b8c6b786844b183b16ffc249c7c5ec704", + "legacy_datatree.zip": "sha256:820cd252047dbf35fa5fb04a9aafee7f7659e0fe4f7d421d69901c57deb6c9d5", + } + + EP = pooch.create( + path=cache_dir, + base_url=base, + version=ver, + registry=registry, + retry_if_failed=1, + ) + + def _unpack(fname, action, pooch_instance): + z = Path(fname) + out = z.parent / z.stem + + print( + "\n[echopype-ci] UNPACK\n" + f" zip file = {z}\n" + f" action = {action}\n" + f" extract_to = {out}\n", + flush=True, + ) + + if action in ("update", "download") or not out.exists(): + from zipfile import ZipFile + + with ZipFile(z, "r") as f: + f.extractall(out) + + # flatten single nested dir if needed + try: + entries = list(out.iterdir()) + if len(entries) == 1 and entries[0].is_dir(): + inner = entries[0] + for child in inner.iterdir(): + target = out / child.name + if not target.exists(): + child.rename(target) + try: + inner.rmdir() + except Exception: + pass + except Exception: + pass + + # Print tree after extraction/flatten + try: + if out.exists(): + print("[echopype-ci] extracted tree (depth ≤ 2):") + for p in sorted(out.glob("*")): + print(f" - {p.name}") + if p.is_dir(): + for q in sorted(p.glob("*")): + print(f" • {q.name}") + except Exception: + pass + + return str(out) + + + for b in bundles: + url = base.format(version=ver) + b + print(f"[echopype-ci] fetching bundle: {b}") + print(f"[echopype-ci] → URL: {url}") + EP.fetch(b, processor=_unpack, progressbar=False) + + TEST_DATA_FOLDER = Path(cache_dir) / ver + + print( + "\n[echopype-ci] TEST_DATA_FOLDER\n" + f" path = {TEST_DATA_FOLDER}\n" + f" exists = {TEST_DATA_FOLDER.exists()}\n", + flush=True, + ) + + if TEST_DATA_FOLDER.exists(): + print("[echopype-ci] top-level contents:") + for p in sorted(TEST_DATA_FOLDER.iterdir()): + print(f" - {p.name}") @pytest.fixture(scope="session") @@ -17,6 +147,7 @@ def test_path(): "ROOT": TEST_DATA_FOLDER, "EA640": TEST_DATA_FOLDER / "ea640", "EK60": TEST_DATA_FOLDER / "ek60", + "EK60_CAL_CHUNKS": TEST_DATA_FOLDER / "ek60_calibrate_chunks", "EK60_MISSING_CHANNEL_POWER": TEST_DATA_FOLDER / "ek60_missing_channel_power", "EK80": TEST_DATA_FOLDER / "ek80", "EK80_NEW": TEST_DATA_FOLDER / "ek80_new", @@ -26,6 +157,11 @@ def test_path(): "AZFP": TEST_DATA_FOLDER / "azfp", "AZFP6": TEST_DATA_FOLDER / "azfp6", "AD2CP": TEST_DATA_FOLDER / "ad2cp", + "EK80_MULTIPLEX": TEST_DATA_FOLDER / "ek80_bb_complex_multiplex", + "EK80_DUPE_PING": TEST_DATA_FOLDER / "ek80_duplicate_ping_times", + "EK80_MISSING_SOUND": TEST_DATA_FOLDER / "ek80_missing_sound_velocity_profile", + "EK80_INVALID_ENV": TEST_DATA_FOLDER / "ek80_invalid_env_datagrams", + "EK80_SEQUENCE": TEST_DATA_FOLDER / "ek80_sequence", "EK80_CAL": TEST_DATA_FOLDER / "ek80_bb_with_calibration", "EK80_EXT": TEST_DATA_FOLDER / "ek80_ext", "ECS": TEST_DATA_FOLDER / "ecs", @@ -40,3 +176,56 @@ def minio_bucket(): key="minioadmin", secret="minioadmin", ) + +# recap of the errors, and failures +def pytest_terminal_summary(terminalreporter, exitstatus, config): + tr = terminalreporter + # Counts + sections = [ + ("passed", "green"), + ("xfailed", "yellow"), + ("xpassed", "yellow"), + ("skipped", "yellow"), + ("failed", "red"), + ("error", "red"), + ] + + tr.write_line("\n─ Test Outcome Summary ─", bold=True) + for key, color in sections: + reports = tr.getreports(key) or [] + if reports: + tr.write_line(f"{key:>8}: {len(reports)}", + **{color: True}) + + # Warnings count (pytest already prints a warnings summary block) + w_reports = tr.stats.get("warnings", []) + if w_reports: + tr.write_line(f"{'warnings':>8}: {len(w_reports)}", yellow=True) + + # List skipped with reasons (short) + skipped = tr.stats.get("skipped", []) + if skipped: + tr.write_line("\nSkipped tests:", yellow=True, bold=True) + for rep in skipped[:20]: + reason = "" + try: + # Best effort to extract reason text + reason = getattr(rep, "longrepr", "") or "" + reason = getattr(reason, "reprcrash", None) and reason.reprcrash.message or str(reason) + except Exception: + pass + tr.write_line(f" • {rep.nodeid}" + (f" — {reason}" if reason else "")) + + if len(skipped) > 20: + tr.write_line(f" … and {len(skipped) - 20} more", yellow=True) + + # List xfailed (expected fails) + xfailed = tr.stats.get("xfailed", []) + if xfailed: + tr.write_line("\nExpected failures (xfail):", yellow=True, bold=True) + for rep in xfailed[:20]: + tr.write_line(f" • {rep.nodeid}") + if len(xfailed) > 20: + tr.write_line(f" … and {len(xfailed) - 20} more", yellow=True) + + tr.write_line("") # trailing newline diff --git a/echopype/tests/consolidate/test_add_depth.py b/echopype/tests/consolidate/test_add_depth.py index 3d75f7fcd..b12e6a5ee 100644 --- a/echopype/tests/consolidate/test_add_depth.py +++ b/echopype/tests/consolidate/test_add_depth.py @@ -3,6 +3,7 @@ import pandas as pd import numpy as np from scipy.spatial.transform import Rotation as R +import os from echopype.utils.align import align_to_ping_time @@ -11,6 +12,18 @@ ek_use_platform_vertical_offsets, ek_use_platform_angles, ek_use_beam_angles ) +# ---- model-root fixtures +@pytest.fixture(scope="module") +def ek60_path(test_path): + return test_path["EK60"] + +@pytest.fixture(scope="module") +def ek80_path(test_path): + return test_path["EK80"] + +@pytest.fixture(scope="module") +def azfp_path(test_path): + return test_path["AZFP"] def _build_ds_Sv(channel, range_sample, ping_time, sample_interval): return xr.Dataset( @@ -175,35 +188,27 @@ def test_warning_zero_vector(caplog): @pytest.mark.integration -@pytest.mark.parametrize("file, sonar_model, compute_Sv_kwargs", [ - ( - "echopype/test_data/ek60/NBP_B050N-D20180118-T090228.raw", - "EK60", - {} - ), - ( - "echopype/test_data/ek60/ncei-wcsd/Summer2017-D20170620-T021537.raw", - "EK60", - {} - ), - ( - "echopype/test_data/ek80/ncei-wcsd/SH1707/Reduced_D20170826-T205615.raw", - "EK80", - {"waveform_mode":"BB", "encode_mode":"complex"} - ), - ( - "echopype/test_data/ek80/ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", - "EK80", - {"waveform_mode":"CW", "encode_mode":"power"} - ) -]) -def test_ek_depth_utils_dims(file, sonar_model, compute_Sv_kwargs): +@pytest.mark.parametrize( + "relpath, sonar_model, compute_Sv_kwargs", + [ + ("NBP_B050N-D20180118-T090228.raw", "EK60", {}), + ("ncei-wcsd/Summer2017-D20170620-T021537.raw", "EK60", {}), + ("ncei-wcsd/SH1707/Reduced_D20170826-T205615.raw", "EK80", {"waveform_mode":"BB", "encode_mode":"complex"}), + ("ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", "EK80", {"waveform_mode":"CW", "encode_mode":"power"}), + ], +) + +def test_ek_depth_utils_dims(relpath, sonar_model, compute_Sv_kwargs, ek60_path, ek80_path): """ Tests `ek_use_platform_vertical_offsets`, `ek_use_platform_angles`, and `ek_use_beam_angles` for correct dimensions. """ - # Open EK Raw file and Compute Sv - ed = ep.open_raw(file, sonar_model=sonar_model) + base = ek60_path if sonar_model == "EK60" else ek80_path + raw_file = base / relpath + if not os.path.isfile(raw_file): + pytest.skip(f"Missing test RAW: {raw_file}") + + ed = ep.open_raw(raw_file, sonar_model=sonar_model) ds_Sv = ep.calibrate.compute_Sv(ed, **compute_Sv_kwargs) # Check dimensions for using EK platform vertical offsets to compute @@ -228,18 +233,19 @@ def test_ek_depth_utils_dims(file, sonar_model, compute_Sv_kwargs): @pytest.mark.integration -def test_ek_depth_utils_group_variable_NaNs_logger_warnings(caplog): +def test_ek_depth_utils_group_variable_NaNs_logger_warnings(caplog, ek80_path): """ Tests `ek_use_platform_vertical_offsets`, `ek_use_platform_angles`, and `ek_use_beam_angles` for correct logger warnings when NaNs exist in group variables. """ + raw_file = ek80_path / "ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw" + if not os.path.isfile(raw_file): + pytest.skip(f"Missing EK80 RAW: {raw_file}") + # Open EK Raw file and Compute Sv - ed = ep.open_raw( - "echopype/test_data/ek80/ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", - sonar_model="EK80" - ) - ds_Sv = ep.calibrate.compute_Sv(ed, **{"waveform_mode":"CW", "encode_mode":"power"}) + ed = ep.open_raw(raw_file, sonar_model="EK80") + ds_Sv = ep.calibrate.compute_Sv(ed, waveform_mode="CW", encode_mode="power") # Set first index of group variables to NaN ed["Platform"]["water_level"].values = np.nan # Is a scalar @@ -285,17 +291,18 @@ def test_ek_depth_utils_group_variable_NaNs_logger_warnings(caplog): @pytest.mark.integration -def test_add_depth_tilt_depth_use_arg_logger_warnings(caplog): +def test_add_depth_tilt_depth_use_arg_logger_warnings(caplog, ek80_path): """ Tests warnings when `tilt` and `depth_offset` are being passed in when other `use_*` arguments are passed in as `True`. """ + raw_file = ek80_path / "ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw" + if not os.path.isfile(raw_file): + pytest.skip(f"Missing EK80 RAW: {raw_file}") + # Open EK Raw file and Compute Sv - ed = ep.open_raw( - "echopype/test_data/ek80/ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", - sonar_model="EK80" - ) - ds_Sv = ep.calibrate.compute_Sv(ed, **{"waveform_mode":"CW", "encode_mode":"power"}) + ed = ep.open_raw(raw_file, sonar_model="EK80") + ds_Sv = ep.calibrate.compute_Sv(ed, waveform_mode="CW", encode_mode="power") # Turn on logger verbosity ep.utils.log.verbose(override=False) @@ -359,14 +366,15 @@ def test_add_depth_without_echodata(): @pytest.mark.integration -def test_add_depth_errors(): +def test_add_depth_errors(ek80_path): """Check if all `add_depth` errors are raised appropriately.""" - # Open EK80 Raw file and Compute Sv - ed = ep.open_raw( - "echopype/test_data/ek80/ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", - sonar_model="EK80" - ) - ds_Sv = ep.calibrate.compute_Sv(ed, **{"waveform_mode":"CW", "encode_mode":"power"}) + raw_file = ek80_path / "ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw" + if not os.path.isfile(raw_file): + pytest.skip(f"Missing EK80 RAW: {raw_file}") + + # Open EK Raw file and Compute Sv + ed = ep.open_raw(raw_file, sonar_model="EK80") + ds_Sv = ep.calibrate.compute_Sv(ed, waveform_mode="CW", encode_mode="power") # Test that all three errors are called: with pytest.raises(ValueError, match=( @@ -388,32 +396,24 @@ def test_add_depth_errors(): @pytest.mark.integration -@pytest.mark.parametrize("file, sonar_model, compute_Sv_kwargs", [ - ( - "echopype/test_data/ek60/NBP_B050N-D20180118-T090228.raw", - "EK60", - {} - ), - ( - "echopype/test_data/ek60/ncei-wcsd/Summer2017-D20170620-T021537.raw", - "EK60", - {} - ), - ( - "echopype/test_data/ek80/ncei-wcsd/SH1707/Reduced_D20170826-T205615.raw", - "EK80", - {"waveform_mode":"BB", "encode_mode":"complex"} - ), - ( - "echopype/test_data/ek80/ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", - "EK80", - {"waveform_mode":"CW", "encode_mode":"power"} - ) -]) -def test_add_depth_EK_with_platform_vertical_offsets(file, sonar_model, compute_Sv_kwargs): +@pytest.mark.parametrize( + "relpath, sonar_model, compute_Sv_kwargs", + [ + ("NBP_B050N-D20180118-T090228.raw", "EK60", {}), + ("ncei-wcsd/Summer2017-D20170620-T021537.raw", "EK60", {}), + ("ncei-wcsd/SH1707/Reduced_D20170826-T205615.raw", "EK80", {"waveform_mode":"BB", "encode_mode":"complex"}), + ("ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", "EK80", {"waveform_mode":"CW", "encode_mode":"power"}), + ], +) +def test_add_depth_EK_with_platform_vertical_offsets(relpath, sonar_model, compute_Sv_kwargs, ek60_path, ek80_path): """Test `depth` values when using EK Platform vertical offset values to compute it.""" - # Open EK Raw file and Compute Sv - ed = ep.open_raw(file, sonar_model=sonar_model) + + base = ek60_path if sonar_model == "EK60" else ek80_path + raw_file = base / relpath + if not os.path.isfile(raw_file): + pytest.skip(f"Missing test RAW: {raw_file}") + + ed = ep.open_raw(raw_file, sonar_model=sonar_model) ds_Sv = ep.calibrate.compute_Sv(ed, **compute_Sv_kwargs) # Subset ds_Sv to include only first 5 `range_sample` coordinates @@ -448,32 +448,21 @@ def test_add_depth_EK_with_platform_vertical_offsets(file, sonar_model, compute_ @pytest.mark.integration -@pytest.mark.parametrize("file, sonar_model, compute_Sv_kwargs", [ - ( - "echopype/test_data/ek60/NBP_B050N-D20180118-T090228.raw", - "EK60", - {} - ), - ( - "echopype/test_data/ek60/ncei-wcsd/Summer2017-D20170620-T021537.raw", - "EK60", - {} - ), - ( - "echopype/test_data/ek80/ncei-wcsd/SH1707/Reduced_D20170826-T205615.raw", - "EK80", - {"waveform_mode":"BB", "encode_mode":"complex"} - ), - ( - "echopype/test_data/ek80/ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", - "EK80", - {"waveform_mode":"CW", "encode_mode":"power"} - ) +@pytest.mark.parametrize("subpath, sonar_model, compute_Sv_kwargs", [ + ("NBP_B050N-D20180118-T090228.raw", "EK60", {}), + ("ncei-wcsd/Summer2017-D20170620-T021537.raw", "EK60", {}), + ("ncei-wcsd/SH1707/Reduced_D20170826-T205615.raw", "EK80", {"waveform_mode": "BB", "encode_mode": "complex"}), + ("ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", "EK80", {"waveform_mode": "CW", "encode_mode": "power"}), ]) -def test_add_depth_EK_with_platform_angles(file, sonar_model, compute_Sv_kwargs): +def test_add_depth_EK_with_platform_angles(subpath, sonar_model, compute_Sv_kwargs, ek60_path, ek80_path): """Test `depth` values when using EK Platform angles to compute it.""" + base = ek60_path if sonar_model == "EK60" else ek80_path + raw_file = base / subpath + if not os.path.isfile(raw_file): + pytest.skip(f"Missing RAW file for {sonar_model}: {raw_file}") + # Open EK Raw file and Compute Sv - ed = ep.open_raw(file, sonar_model=sonar_model) + ed = ep.open_raw(raw_file, sonar_model=sonar_model) ds_Sv = ep.calibrate.compute_Sv(ed, **compute_Sv_kwargs) # Replace any Beam Angle NaN values with 0 @@ -502,32 +491,21 @@ def test_add_depth_EK_with_platform_angles(file, sonar_model, compute_Sv_kwargs) @pytest.mark.integration -@pytest.mark.parametrize("file, sonar_model, compute_Sv_kwargs", [ - ( - "echopype/test_data/ek60/NBP_B050N-D20180118-T090228.raw", - "EK60", - {} - ), - ( - "echopype/test_data/ek60/ncei-wcsd/Summer2017-D20170620-T021537.raw", - "EK60", - {} - ), - ( - "echopype/test_data/ek80/ncei-wcsd/SH1707/Reduced_D20170826-T205615.raw", - "EK80", - {"waveform_mode":"BB", "encode_mode":"complex"} - ), - ( - "echopype/test_data/ek80/ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", - "EK80", - {"waveform_mode":"CW", "encode_mode":"power"} - ) +@pytest.mark.parametrize("subpath, sonar_model, compute_Sv_kwargs", [ + ("NBP_B050N-D20180118-T090228.raw", "EK60", {}), + ("ncei-wcsd/Summer2017-D20170620-T021537.raw", "EK60", {}), + ("ncei-wcsd/SH1707/Reduced_D20170826-T205615.raw", "EK80", {"waveform_mode": "BB", "encode_mode": "complex"}), + ("ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", "EK80", {"waveform_mode": "CW", "encode_mode": "power"}), ]) -def test_add_depth_EK_with_beam_angles(file, sonar_model, compute_Sv_kwargs): +def test_add_depth_EK_with_beam_angles(subpath, sonar_model, compute_Sv_kwargs, ek60_path, ek80_path): """Test `depth` values when using EK Beam angles to compute it.""" + base = ek60_path if sonar_model == "EK60" else ek80_path + raw_file = base / subpath + if not os.path.isfile(raw_file): + pytest.skip(f"Missing RAW file for {sonar_model}: {raw_file}") + # Open EK Raw file and Compute Sv - ed = ep.open_raw(file, sonar_model=sonar_model) + ed = ep.open_raw(raw_file, sonar_model=sonar_model) ds_Sv = ep.calibrate.compute_Sv(ed, **compute_Sv_kwargs) # Replace Beam Angle NaN values @@ -548,6 +526,31 @@ def test_add_depth_EK_with_beam_angles(file, sonar_model, compute_Sv_kwargs): # Compute echo range scaling values echo_range_scaling = ek_use_beam_angles(ed["Sonar/Beam_group1"]) + # Check if depth is equal to echo range scaling value * echo range + assert np.allclose( + ds_Sv_with_depth["depth"].data, + (echo_range_scaling * ds_Sv["echo_range"]).transpose("channel", "ping_time", "range_sample").data, + equal_nan=True, + ) + + # Replace Beam Angle NaN values + ed["Sonar/Beam_group1"]["beam_direction_x"].values = ed["Sonar/Beam_group1"]["beam_direction_x"].fillna(0).values + ed["Sonar/Beam_group1"]["beam_direction_y"].values = ed["Sonar/Beam_group1"]["beam_direction_y"].fillna(0).values + ed["Sonar/Beam_group1"]["beam_direction_z"].values = ed["Sonar/Beam_group1"]["beam_direction_z"].fillna(1).values + + # Compute `depth` using beam angle values + ds_Sv_with_depth = ep.consolidate.add_depth(ds_Sv, ed, use_beam_angles=True) + + # Check history attribute + history_attribute = ds_Sv_with_depth["depth"].attrs["history"] + history_attribute_without_time = history_attribute[32:] + assert history_attribute_without_time == ( + ". `depth` calculated using: Sv `echo_range`, Echodata `Beam_group1` Angles." + ) + + # Compute echo range scaling values + echo_range_scaling = ek_use_beam_angles(ed["Sonar/Beam_group1"]) + # Check if depth is equal to echo range scaling value * echo range assert np.allclose( ds_Sv_with_depth["depth"].data, @@ -555,33 +558,72 @@ def test_add_depth_EK_with_beam_angles(file, sonar_model, compute_Sv_kwargs): equal_nan=True ) - @pytest.mark.integration -@pytest.mark.parametrize("file, sonar_model, compute_Sv_kwargs, expected_beam_group_name", [ - ( - "echopype/test_data/ek80/Summer2018--D20180905-T033113.raw", - "EK80", - {"waveform_mode":"BB", "encode_mode":"complex"}, - "Beam_group1" - ), - ( - "echopype/test_data/ek80/Summer2018--D20180905-T033113.raw", - "EK80", - {"waveform_mode":"CW", "encode_mode":"power"}, - "Beam_group2" - ) +@pytest.mark.parametrize("file, sonar_model, compute_Sv_kwargs", [ + ("NBP_B050N-D20180118-T090228.raw", "EK60", {}), + ("ncei-wcsd/SH1707/Reduced_D20170826-T205615.raw", "EK80", {"waveform_mode": "BB", "encode_mode": "complex"}), + ("ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", "EK80", {"waveform_mode": "CW", "encode_mode": "power"}) ]) +def test_add_depth_with_dim_swap(file, sonar_model, compute_Sv_kwargs, ek80_path, ek60_path): + """ + Test adding depth to Sv dataset after swapping dimension/coordinate + from channel to frequency_nominal. + Asserts that the output dataset has swapped channel dim to frequency_nominal + and contains the depth variable. + """ + if sonar_model == "EK60": + ed = ep.open_raw(ek60_path / file, sonar_model=sonar_model) + else: + ed = ep.open_raw(ek80_path / file, sonar_model=sonar_model) + + ds_Sv = ep.calibrate.compute_Sv(ed, **compute_Sv_kwargs) + + ds_Sv = ep.consolidate.swap_dims_channel_frequency(ds_Sv) + + # swap dims in beam_group to test with dim_0 = frequency_nominal + ed["Sonar/Beam_group1"] = ep.consolidate.swap_dims_channel_frequency(ed["Sonar/Beam_group1"]) + + # Replace Beam Angle NaN values + ed["Sonar/Beam_group1"]["beam_direction_x"].values = ed["Sonar/Beam_group1"]["beam_direction_x"].fillna(0).values + ed["Sonar/Beam_group1"]["beam_direction_y"].values = ed["Sonar/Beam_group1"]["beam_direction_y"].fillna( 0).values + ed["Sonar/Beam_group1"]["beam_direction_z"].values = ed["Sonar/Beam_group1"]["beam_direction_z"].fillna(1).values + + ds_Sv_with_depth = ep.consolidate.add_depth(ds_Sv, ed, use_beam_angles=True) + + # Check that channel dim has been swapped to frequency_nominal + assert "channel" not in ds_Sv_with_depth.sizes + assert "frequency_nominal" in ds_Sv_with_depth.sizes + # Check that depth has been added + assert "depth" in ds_Sv_with_depth.data_vars + + +@pytest.mark.integration +@pytest.mark.parametrize( + "subpath, compute_Sv_kwargs, expected_beam_group_name", + [ + ("Summer2018--D20180905-T033113.raw", + {"waveform_mode": "BB", "encode_mode": "complex"}, + "Beam_group1"), + ("Summer2018--D20180905-T033113.raw", + {"waveform_mode": "CW", "encode_mode": "power"}, + "Beam_group2"), + ], +) def test_add_depth_EK_with_beam_angles_with_different_beam_groups( - file, sonar_model, compute_Sv_kwargs, expected_beam_group_name + subpath, compute_Sv_kwargs, expected_beam_group_name, ek80_path ): """ Test `depth` channel when using EK Beam angles from two separate calibrated - Sv datasets (that are from the same raw file) using two differing pairs of - calibration key word arguments. The two tests should correspond to different - beam groups i.e. beam group 1 and beam group 2. + Sv datasets (same raw file) with different calibration kwargs, yielding + different beam groups. """ + sonar_model = "EK80" + raw_file = ek80_path / subpath + if not os.path.isfile(raw_file): + pytest.skip(f"Missing EK80 RAW: {raw_file}") + # Open EK Raw file and Compute Sv - ed = ep.open_raw(file, sonar_model=sonar_model) + ed = ep.open_raw(raw_file, sonar_model=sonar_model) ds_Sv = ep.calibrate.compute_Sv(ed, **compute_Sv_kwargs) # Compute `depth` using beam angle values @@ -594,24 +636,30 @@ def test_add_depth_EK_with_beam_angles_with_different_beam_groups( f". `depth` calculated using: Sv `echo_range`, Echodata `{expected_beam_group_name}` Angles." ) - @pytest.mark.integration -def test_add_depth_with_external_glider_depth_and_tilt_array(): +def test_add_depth_with_external_glider_depth_and_tilt_array(azfp_path): """ Test add_depth with external glider depth offset and tilt array data. """ + # Define paths + raw_file = azfp_path / "rutgers_glider_external_nc/18011107.01A" + xml_file = azfp_path / "rutgers_glider_external_nc/18011107.XML" + glider_file = azfp_path / "rutgers_glider_external_nc/ru32-20180109T0531-profile-sci-delayed-subset.nc" + + # Skip if missing + for p in [raw_file, xml_file, glider_file]: + if not os.path.isfile(p): + pytest.skip(f"Missing required AZFP test data file: {p}") + # Open RAW ed = ep.open_raw( - raw_file="echopype/test_data/azfp/rutgers_glider_external_nc/18011107.01A", - xml_path="echopype/test_data/azfp/rutgers_glider_external_nc/18011107.XML", - sonar_model="azfp" + raw_file=raw_file, + xml_path=xml_file, + sonar_model="azfp", ) # Open external glider dataset - glider_ds = xr.open_dataset( - "echopype/test_data/azfp/rutgers_glider_external_nc/ru32-20180109T0531-profile-sci-delayed-subset.nc", - engine="netcdf4" - ) + glider_ds = xr.open_dataset(glider_file, engine="netcdf4") # Grab external environment parameters env_params_means = {} @@ -663,18 +711,20 @@ def test_add_depth_with_external_glider_depth_and_tilt_array(): @pytest.mark.unit -def test_multi_dim_depth_offset_and_tilt_array_error(): +def test_multi_dim_depth_offset_and_tilt_array_error(ek80_path): """ Test that the correct `ValueError`s are raised when a multi-dimensional array is passed into `add_depth` for the `depth_offset` and `tilt` arguments. """ + # Define test file path + raw_file = ek80_path / "ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw" + if not os.path.isfile(raw_file): + pytest.skip(f"Missing EK80 RAW: {raw_file}") + # Open EK Raw file and Compute Sv - ed = ep.open_raw( - "echopype/test_data/ek80/ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131621.raw", - sonar_model="EK80" - ) - ds_Sv = ep.calibrate.compute_Sv(ed, **{"waveform_mode":"CW", "encode_mode":"power"}) + ed = ep.open_raw(raw_file, sonar_model="EK80") + ds_Sv = ep.calibrate.compute_Sv(ed, waveform_mode="CW", encode_mode="power") # Multi-dimensional mock array multi_dim_da = xr.DataArray( diff --git a/echopype/tests/consolidate/test_add_location.py b/echopype/tests/consolidate/test_add_location.py index a22fe414b..05159465a 100644 --- a/echopype/tests/consolidate/test_add_location.py +++ b/echopype/tests/consolidate/test_add_location.py @@ -194,6 +194,73 @@ def _tests(ds_test, location_type, nmea_sentence=None): ds_sel = ep.consolidate.add_location(ds=ds, echodata=ed, nmea_sentence="GGA") _tests(ds_sel, location_type, nmea_sentence="GGA") +@pytest.mark.integration +@pytest.mark.parametrize( + ["sonar_model", "path_model", "raw_and_xml_paths", "lat_lon_name_dict", "extras"], + [ + ( + "AZFP", + "AZFP", + ("17082117.01A", "17041823.XML"), + {"lat_name": "latitude", "lon_name": "longitude"}, + {'longitude': -60.0, 'latitude': 45.0, 'salinity': 27.9, 'pressure': 59}, + ), + ], +) +def test_add_location_with_dim_swap( + sonar_model, + path_model, + raw_and_xml_paths, + lat_lon_name_dict, + extras, + test_path +): + """ + Test adding location to Sv dataset after swapping dimension/coordinate + from channel to frequency_nominal. + Asserts that the output dataset has swapped channel dim to frequency_nominal + and contains the latitude and longitude variable. + """ + + raw_path = test_path[path_model] / raw_and_xml_paths[0] + xml_path = test_path[path_model] / raw_and_xml_paths[1] + + ed = ep.open_raw(raw_path, xml_path=xml_path, sonar_model=sonar_model) + point_ds = xr.Dataset( + { + lat_lon_name_dict["lat_name"]: (["time"], np.array([float(extras['latitude'])])), + lat_lon_name_dict["lon_name"]: (["time"], np.array([float(extras['longitude'])])), + }, + coords={ + "time": (["time"], np.array([ed["Sonar/Beam_group1"]["ping_time"].values.min()])) + }, + ) + ed.update_platform( + point_ds, + variable_mappings={ + lat_lon_name_dict["lat_name"]: lat_lon_name_dict["lat_name"], + lat_lon_name_dict["lon_name"]: lat_lon_name_dict["lon_name"] + } + ) + + env_params = { + "temperature": ed["Environment"]["temperature"].values.mean(), + "salinity": extras["salinity"], + "pressure": extras["pressure"], + } + + ds = ep.calibrate.compute_Sv(echodata=ed, env_params=env_params) + + ds = ep.consolidate.swap_dims_channel_frequency(ds) + + ds_all = ep.consolidate.add_location(ds=ds, echodata=ed) + + # Check that channel dim has been swapped to frequency_nominal + assert "channel" not in ds_all.sizes + assert "frequency_nominal" in ds_all.sizes + # Check that latitude and longitude have been added + assert "latitude" in ds_all.data_vars + assert "longitude" in ds_all.data_vars @pytest.mark.integration @pytest.mark.parametrize( diff --git a/echopype/tests/consolidate/test_consolidate_integration.py b/echopype/tests/consolidate/test_consolidate_integration.py index 03341e07d..a17bb3ca5 100644 --- a/echopype/tests/consolidate/test_consolidate_integration.py +++ b/echopype/tests/consolidate/test_consolidate_integration.py @@ -286,6 +286,68 @@ def test_add_splitbeam_angle(sonar_model, test_path_key, raw_file_name, test_pat # remove the temporary directory, if it was created temp_dir.cleanup() +@pytest.mark.parametrize( + ("sonar_model", "test_path_key", "raw_file_name", "paths_to_echoview_mat"), + [ + # ek60_CW_power + ( + "EK60", "EK60", "DY1801_EK60-D20180211-T164025.raw", + [ + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T1.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T2.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T3.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T4.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T5.mat' + ], + ), + # ek80_CW_power + ( + "EK80", "EK80", "Summer2018--D20180905-T033113.raw", + [ + 'splitbeam/Summer2018--D20180905-T033113_angles_T2.mat', + 'splitbeam/Summer2018--D20180905-T033113_angles_T1.mat', + ], + ), + ], + ids=[ + "ek60_CW_power", + "ek80_CW_power", + ], +) +def test_add_splitbeam_angle_with_dim_swap(sonar_model, test_path_key, raw_file_name, test_path, paths_to_echoview_mat): + """ + Test adding split-beam angle to Sv dataset after swapping dimension/coordinate + from channel to frequency_nominal. + Asserts that the output dataset has swapped channel dim to frequency_nominal + and contains the split-beam angle variables. + """ + + ed = ep.open_raw(test_path[test_path_key] / raw_file_name, sonar_model=sonar_model) + + waveform_mode = "CW" + encode_mode = "power" + + ds_Sv = ep.calibrate.compute_Sv(ed, waveform_mode=waveform_mode, encode_mode=encode_mode) + + ds_Sv = ep.consolidate.swap_dims_channel_frequency(ds_Sv) + + # swap dims in beam_groups to test with dim_0 = frequency_nominal + ed["Sonar/Beam_group1"] = ep.consolidate.swap_dims_channel_frequency(ed["Sonar/Beam_group1"]) + if ed["Sonar"].sizes["beam_group"] > 1: + ed["Sonar/Beam_group2"] = ep.consolidate.swap_dims_channel_frequency(ed["Sonar/Beam_group2"]) + + ds_Sv = ep.consolidate.add_splitbeam_angle(source_Sv=ds_Sv, echodata=ed, + waveform_mode=waveform_mode, + encode_mode=encode_mode, + to_disk=False) + + # Check that channel dim has been swapped to frequency_nominal + assert "channel" not in ds_Sv.sizes + assert "frequency_nominal" in ds_Sv.sizes + # Check that split-beam angles were added to the dataset + assert "angle_alongship" in ds_Sv.data_vars + assert "angle_athwartship" in ds_Sv.data_vars + def test_add_splitbeam_angle_BB_pc(test_path): diff --git a/echopype/tests/convert/test_convert_ad2cp.py b/echopype/tests/convert/test_convert_ad2cp.py index 07386bdad..8addaa2fe 100644 --- a/echopype/tests/convert/test_convert_ad2cp.py +++ b/echopype/tests/convert/test_convert_ad2cp.py @@ -14,14 +14,16 @@ from pathlib import Path from echopype import open_raw, open_converted -from echopype.testing import TEST_DATA_FOLDER +@pytest.fixture(scope="session") +def ad2cp_path(test_path): + return test_path["AD2CP"] + @pytest.fixture def ocean_contour_export_dir(test_path): return test_path["AD2CP"] / "ocean-contour" - @pytest.fixture def ocean_contour_export_076_dir(ocean_contour_export_dir): return ocean_contour_export_dir / "076" @@ -38,13 +40,17 @@ def output_dir(): def pytest_generate_tests(metafunc): - ad2cp_path = TEST_DATA_FOLDER / "ad2cp" - test_file_dir = ( - ad2cp_path / "normal" - ) # "normal" files do not have IQ samples - raw_test_file_dir = ad2cp_path / "raw" # "raw" files contain IQ samples - ad2cp_files = test_file_dir.glob("**/*.ad2cp") - raw_ad2cp_files = raw_test_file_dir.glob("**/*.ad2cp") + # Import the module echopype.tests + # not the fixture as they are not meant to be called directly, + from echopype.tests import conftest as ct + + ad2cp_root = ct.TEST_DATA_FOLDER / "ad2cp" + test_file_dir = ad2cp_root / "normal" # no IQ samples + raw_test_file_dir = ad2cp_root / "raw" # has IQ samples + + ad2cp_files = sorted(test_file_dir.glob("**/*.ad2cp")) + raw_ad2cp_files = sorted(raw_test_file_dir.glob("**/*.ad2cp")) + if "filepath" in metafunc.fixturenames: metafunc.parametrize( argnames="filepath", @@ -75,12 +81,11 @@ def absolute_tolerance(): return 1e-6 -def test_convert(filepath, output_dir): - with TemporaryDirectory() as tmpdir: - output_dir = Path(tmpdir + output_dir) - print("converting", filepath) - echodata = open_raw(raw_file=str(filepath), sonar_model="AD2CP") - echodata.to_netcdf(save_path=output_dir) +def test_convert(filepath, output_dir, tmp_path): + output_dir = tmp_path / output_dir.strip("/") + print("converting", filepath) + echodata = open_raw(raw_file=str(filepath), sonar_model="AD2CP") + echodata.to_netcdf(save_path=output_dir) def test_convert_raw( @@ -89,22 +94,21 @@ def test_convert_raw( ocean_contour_export_090_dir, ocean_contour_export_076_dir, absolute_tolerance, + tmp_path, ): - with TemporaryDirectory() as tmpdir: - output_dir = Path(tmpdir + output_dir) - - print("converting raw", filepath_raw) - echodata = open_raw(raw_file=str(filepath_raw), sonar_model="AD2CP") - echodata.to_netcdf(save_path=output_dir) - - _check_raw_output( - filepath_raw, - output_dir, - ocean_contour_export_090_dir, - ocean_contour_export_076_dir, - absolute_tolerance, - ) + output_dir = tmp_path / output_dir.strip("/") + + print("converting raw", filepath_raw) + echodata = open_raw(raw_file=str(filepath_raw), sonar_model="AD2CP") + echodata.to_netcdf(save_path=output_dir) + _check_raw_output( + filepath_raw, + output_dir, + ocean_contour_export_090_dir, + ocean_contour_export_076_dir, + absolute_tolerance, + ) def _check_raw_output( filepath_raw, @@ -261,4 +265,4 @@ def _check_raw_output( base["Data_Q"].data.T.flatten(), atol=absolute_tolerance, ) - base.close() + base.close() \ No newline at end of file diff --git a/echopype/tests/convert/test_convert_ek.py b/echopype/tests/convert/test_convert_ek.py index b73163021..42a45c383 100644 --- a/echopype/tests/convert/test_convert_ek.py +++ b/echopype/tests/convert/test_convert_ek.py @@ -6,10 +6,17 @@ from echopype.convert.utils.ek_raw_io import RawSimradFile, SimradEOF from echopype.convert.parse_base import ParseEK +@pytest.fixture +def ek60_path(test_path): + return test_path["EK60"] + +@pytest.fixture +def ek80_path(test_path): + return test_path["EK80"] def expected_array_shape(file, datagram_type, datagram_item): """Extract array shape from user-specified parsed datagram type and item.""" - fid = RawSimradFile(file.replace("raw", datagram_type)) + fid = RawSimradFile(str(file).replace("raw", datagram_type)) fid.read(1) datagram_item_list = [] while True: @@ -21,19 +28,19 @@ def expected_array_shape(file, datagram_type, datagram_item): @pytest.mark.integration -def test_convert_ek60_with_missing_bot_idx_file(): +def test_convert_ek60_with_missing_bot_idx_file(ek60_path, ek80_path): """ Check appropriate FileNotFoundError when attempting to parse .BOT and IDX file that do not exist in the same folder for which the .RAW file exists.""" with pytest.raises(FileNotFoundError): open_raw( - "echopype/test_data/ek60/ncei-wcsd/SH1701/TEST-D20170114-T202932.raw", + ek60_path / "ncei-wcsd/SH1701/TEST-D20170114-T202932.raw", sonar_model="EK60", include_bot=True, ) with pytest.raises(FileNotFoundError): open_raw( - "echopype/test_data/ek80/ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131325.raw", + ek80_path / "ncei-wcsd/SH2106/EK80/Reduced_Hake-D20210701-T131325.raw", sonar_model="EK80", include_idx=True, ) @@ -43,20 +50,18 @@ def test_convert_ek60_with_missing_bot_idx_file(): @pytest.mark.parametrize( "file, sonar_model", [ - ("echopype/test_data/ek60/idx_bot/Summer2017-D20170620-T011027.raw", "EK60"), - ("echopype/test_data/ek60/idx_bot/Summer2017-D20170707-T150923.raw", "EK60"), - ("echopype/test_data/ek80/idx_bot/Hake-D20230711-T181910.raw", "EK80"), - ("echopype/test_data/ek80/idx_bot/Hake-D20230711-T182702.raw", "EK80"), + ("idx_bot/Summer2017-D20170620-T011027.raw", "EK60"), + ("idx_bot/Summer2017-D20170707-T150923.raw", "EK60"), + ("idx_bot/Hake-D20230711-T181910.raw", "EK80"), + ("idx_bot/Hake-D20230711-T182702.raw", "EK80"), ] ) -def test_convert_ek_with_bot_file(file, sonar_model): +def test_convert_ek_with_bot_file(file, sonar_model, ek60_path, ek80_path): """Check variable dimensions, time encodings, and attributes when BOT file is parsed.""" # Open Raw and Parse BOT - ed = open_raw( - file, - sonar_model=sonar_model, - include_bot=True, - ) + path = ek60_path if sonar_model == "EK60" else ek80_path + file = path / file + ed = open_raw(file, sonar_model=sonar_model, include_bot=True) # Check data variable shape seafloor_depth_da = ed["Vendor_specific"]["detected_seafloor_depth"] @@ -89,20 +94,18 @@ def test_convert_ek_with_bot_file(file, sonar_model): @pytest.mark.parametrize( "file, sonar_model", [ - ("echopype/test_data/ek60/idx_bot/Summer2017-D20170620-T011027.raw", "EK60"), - ("echopype/test_data/ek60/idx_bot/Summer2017-D20170707-T150923.raw", "EK60"), - ("echopype/test_data/ek80/idx_bot/Hake-D20230711-T181910.raw", "EK80"), - ("echopype/test_data/ek80/idx_bot/Hake-D20230711-T182702.raw", "EK80"), + ("idx_bot/Summer2017-D20170620-T011027.raw", "EK60"), + ("idx_bot/Summer2017-D20170707-T150923.raw", "EK60"), + ("idx_bot/Hake-D20230711-T181910.raw", "EK80"), + ("idx_bot/Hake-D20230711-T182702.raw", "EK80"), ] ) -def test_convert_ek_with_idx_file(file, sonar_model): +def test_convert_ek_with_idx_file(file, sonar_model, ek60_path, ek80_path): """Check variable dimensions and attributes when IDX file is parsed.""" # Open Raw and Parse IDX - ed = open_raw( - file, - sonar_model=sonar_model, - include_idx=True, - ) + path = ek60_path if sonar_model == "EK60" else ek80_path + file = path / file + ed = open_raw(file, sonar_model=sonar_model, include_idx=True) platform = ed["Platform"] # Check data variable lengths diff --git a/echopype/tests/convert/test_convert_ek60.py b/echopype/tests/convert/test_convert_ek60.py index 054e7fef6..a44f22cfa 100644 --- a/echopype/tests/convert/test_convert_ek60.py +++ b/echopype/tests/convert/test_convert_ek60.py @@ -21,17 +21,6 @@ def ek60_missing_channel_power_path(test_path): def es60_path(test_path): return test_path["ES60"] - -# raw_paths = ['./echopype/test_data/ek60/set1/' + file -# for file in os.listdir('./echopype/test_data/ek60/set1')] # 2 range lengths -# raw_path = ['./echopype/test_data/ek60/set2/' + file -# for file in os.listdir('./echopype/test_data/ek60/set2')] # 3 range lengths -# Other data files -# raw_filename = 'data_zplsc/OceanStarr_2017-D20170725-T004612.raw' # OceanStarr 2 channel EK60 -# raw_filename = '../data/DY1801_EK60-D20180211-T164025.raw' # Dyson 5 channel EK60 -# raw_filename = 'data_zplsc/D20180206-T000625.raw # EK80 - - def test_convert_ek60_matlab_raw(ek60_path): """Compare parsed Beam group data with Matlab outputs.""" ek60_raw_path = str( diff --git a/echopype/tests/convert/test_convert_ek80.py b/echopype/tests/convert/test_convert_ek80.py index ae7d5434d..34d8ca5e2 100644 --- a/echopype/tests/convert/test_convert_ek80.py +++ b/echopype/tests/convert/test_convert_ek80.py @@ -8,7 +8,6 @@ from echopype import open_raw, open_converted from echopype.calibrate import compute_Sv -from echopype.testing import TEST_DATA_FOLDER from echopype.convert.parse_ek80 import ParseEK80 from echopype.convert.set_groups_ek80 import WIDE_BAND_TRANS, PULSE_COMPRESS, FILTER_IMAG, FILTER_REAL, DECIMATION from echopype.utils import log @@ -19,13 +18,39 @@ def ek80_path(test_path): return test_path["EK80"] +@pytest.fixture +def ek80_dupe_ping_path(test_path): + return test_path["EK80_DUPE_PING"] + +@pytest.fixture +def ek80_missing_sound_path(test_path): + return test_path["EK80_MISSING_SOUND"] + +@pytest.fixture +def ek80_invalid_env_path(test_path): + return test_path["EK80_INVALID_ENV"] + +@pytest.fixture +def ek80_sequence_path(test_path): + return test_path["EK80_SEQUENCE"] + +@pytest.fixture +def ek80_new_path(test_path): + return test_path["EK80_NEW"] def pytest_generate_tests(metafunc): - ek80_new_path = TEST_DATA_FOLDER / "ek80_new" - ek80_new_files = ek80_new_path.glob("**/*.raw") + """Dynamically parameterize tests for EK80 .raw files.""" + from echopype.tests import conftest as ct + + ek80_new_root = ct.TEST_DATA_FOLDER / "ek80_new" + ek80_new_files = sorted(ek80_new_root.glob("**/*.raw")) + + ek80_new_files = sorted(ek80_new_root.glob("**/*.raw")) if "ek80_new_file" in metafunc.fixturenames: metafunc.parametrize( - "ek80_new_file", ek80_new_files, ids=lambda f: str(f.name) + "ek80_new_file", + ek80_new_files, + ids=[f.name for f in ek80_new_files], ) @@ -33,14 +58,6 @@ def pytest_generate_tests(metafunc): def ek80_new_file(request): return request.param - -# raw_path_simrad = ['./echopype/test_data/ek80/simrad/EK80_SimradEcho_WC381_Sequential-D20150513-T090935.raw', -# './echopype/test_data/ek80/simrad/EK80_SimradEcho_WC381_Sequential-D20150513-T091004.raw', -# './echopype/test_data/ek80/simrad/EK80_SimradEcho_WC381_Sequential-D20150513-T091034.raw', -# './echopype/test_data/ek80/simrad/EK80_SimradEcho_WC381_Sequential-D20150513-T091105.raw'] -# raw_paths = ['./echopype/test_data/ek80/Summer2018--D20180905-T033113.raw', -# './echopype/test_data/ek80/Summer2018--D20180905-T033258.raw'] # Multiple files (CW and BB) - def check_env_xml(echodata): # check environment xml datagram @@ -491,19 +508,19 @@ def test_parse_mru0_mru1(ek80_path): @pytest.mark.unit -def test_parse_missing_sound_velocity_profile(): +def test_parse_missing_sound_velocity_profile(ek80_missing_sound_path): """ Tests that RAW files that are missing sound velocity profile values can be converted, saved to Zarr, and opened again. """ # Open RAW ed = open_raw( - "echopype/test_data/ek80_missing_sound_velocity_profile/Hake-D20230701-T073658.raw", + ek80_missing_sound_path / "Hake-D20230701-T073658.raw", sonar_model="EK80" ) # Save RAW to Zarr - save_path = "echopype/test_data/ek80_missing_sound_velocity_profile/test_save.zarr" + save_path = ek80_missing_sound_path / "test_save.zarr" ed.to_zarr(save_path,overwrite=True) # Open Converted @@ -518,7 +535,7 @@ def test_parse_missing_sound_velocity_profile(): @pytest.mark.unit -def test_duplicate_ping_times(caplog): +def test_duplicate_ping_times(caplog, ek80_dupe_ping_path): """ Tests that RAW file with duplicate ping times can be parsed and that the correct warning has been raised. """ @@ -526,7 +543,7 @@ def test_duplicate_ping_times(caplog): log.verbose(override=False) # Open RAW - ed = open_raw("echopype/test_data/ek80_duplicate_ping_times/Hake-D20210913-T130612.raw", sonar_model="EK80") + ed = open_raw(ek80_dupe_ping_path / "Hake-D20210913-T130612.raw", sonar_model="EK80") # Check that there are no ping time duplicates in Beam group assert ed["Sonar/Beam_group1"].equals( @@ -542,7 +559,7 @@ def test_duplicate_ping_times(caplog): @pytest.mark.unit -def test_check_unique_ping_time_duplicates(caplog): +def test_check_unique_ping_time_duplicates(caplog, ek80_dupe_ping_path): """ Checks that `check_unique_ping_time_duplicates` raises a warning when the data for duplicate ping times is not unique. """ @@ -553,7 +570,7 @@ def test_check_unique_ping_time_duplicates(caplog): log.verbose(override=False) # Open duplicate ping time beam dataset - ds_data = xr.open_zarr("echopype/test_data/ek80_duplicate_ping_times/duplicate_beam_ds.zarr") + ds_data = xr.open_zarr(ek80_dupe_ping_path / "duplicate_beam_ds.zarr") # Modify a single entry to ensure that there exists duplicate ping times that do not share the same backscatter data ds_data["backscatter_r"][0,0,0] = 0 @@ -574,7 +591,7 @@ def test_check_unique_ping_time_duplicates(caplog): @pytest.mark.unit -def test_parse_ek80_with_invalid_env_datagrams(): +def test_parse_ek80_with_invalid_env_datagrams(ek80_invalid_env_path): """ Tests parsing EK80 RAW file with invalid environment datagrams. Checks that the EchoData object contains the necessary environment variables for calibration. @@ -582,7 +599,7 @@ def test_parse_ek80_with_invalid_env_datagrams(): # Parse RAW ed = open_raw( - "echopype/test_data/ek80_invalid_env_datagrams/SH24-replay-D20240705-T070536.raw", + ek80_invalid_env_path / "SH24-replay-D20240705-T070536.raw", sonar_model="EK80", ) @@ -593,14 +610,14 @@ def test_parse_ek80_with_invalid_env_datagrams(): @pytest.mark.unit -def test_ek80_sonar_all_channel(): +def test_ek80_sonar_all_channel(ek80_new_path): """ Checks that when an EK80 raw file has 2 beam groups, the converted Echodata object contains all channels in the 'channel_all' dimension of the Sonar group. """ # Convert EK80 Raw File ed = open_raw( - raw_file="echopype/test_data/ek80_new/echopype-test-D20211004-T235714.raw", + raw_file= ek80_new_path / "echopype-test-D20211004-T235714.raw", sonar_model="EK80" ) @@ -616,14 +633,14 @@ def test_ek80_sonar_all_channel(): @pytest.mark.unit -def test_ek80_sequence_filter_coeff(): +def test_ek80_sequence_filter_coeff(ek80_sequence_path): """ Checks that filter coefficients are stored properly for EK80 raw files generated with ping sequence. """ # Convert EK80 Raw File ed = open_raw( - raw_file="echopype/test_data/ek80_sequence/three_ensemble-Phase0-D20240506-T053349-0.raw", + raw_file= ek80_sequence_path / "three_ensemble-Phase0-D20240506-T053349-0.raw", sonar_model="EK80" ) diff --git a/echopype/tests/convert/test_convert_source_target_locs.py b/echopype/tests/convert/test_convert_source_target_locs.py index 3f8f6954f..43e1c11dd 100644 --- a/echopype/tests/convert/test_convert_source_target_locs.py +++ b/echopype/tests/convert/test_convert_source_target_locs.py @@ -282,13 +282,32 @@ def test_convert_ek( common_storage_options = minio_bucket output_storage_options = {} input_paths, sonar_model = ek_input_params + ipath = input_paths if isinstance(input_paths, list): + # If the test ever uses multiple HTTP files, wrap each one. + wrapped = [] + for p in input_paths: + if p.startswith(("http://", "https://")): + wrapped.append(f"simplecache::{p}") + else: + wrapped.append(p) + input_paths = wrapped ipath = input_paths[0] - input_storage_options = ( - common_storage_options if ipath.startswith("s3://") else {} - ) + # Handle input sources + if ipath.startswith("s3://"): + input_storage_options = common_storage_options + elif ipath.startswith(("http://", "https://", "simplecache::http://", "simplecache::https://")): + # Make HTTP files seekable across all OSes + if not ipath.startswith("simplecache::"): + ipath = f"simplecache::{ipath}" + # Disable file-change checks to keep it simple/fast in CI + input_storage_options = {} + else: + input_storage_options = {} + + # Handle outputs (S3 only) if output_save_path and output_save_path.startswith("s3://"): output_storage_options = common_storage_options @@ -360,9 +379,15 @@ def test_convert_azfp( common_storage_options = minio_bucket output_storage_options = {} - input_storage_options = ( - common_storage_options if azfp_input_paths.startswith("s3://") else {} - ) + # S3 uses MinIO creds; HTTP must stream to avoid ranged reads on Windows CI + if azfp_input_paths.startswith("s3://"): + input_storage_options = common_storage_options + elif azfp_input_paths.startswith(("http://", "https://")): + azfp_input_paths = f"simplecache::{azfp_input_paths}" + input_storage_options = {} + else: + input_storage_options = {} + if output_save_path and output_save_path.startswith("s3://"): output_storage_options = common_storage_options diff --git a/echopype/tests/echodata/test_echodata.py b/echopype/tests/echodata/test_echodata.py index 01fae4337..40c7dd4f0 100644 --- a/echopype/tests/echodata/test_echodata.py +++ b/echopype/tests/echodata/test_echodata.py @@ -20,6 +20,9 @@ from utils import get_mock_echodata, check_consolidated +@pytest.fixture +def ek60_path(test_path): + return test_path["EK60"] @pytest.fixture(scope="module") def legacy_datatree(test_path): @@ -314,27 +317,21 @@ def test_get_dataset(self, converted_zarr): assert result is None assert isinstance(ed_result, xr.Dataset) - @pytest.mark.parametrize("consolidated", [True, False]) - def test_to_zarr_consolidated(self, mock_echodata, consolidated): + def test_to_zarr_created(self, mock_echodata): """ - Tests to_zarr consolidation. Currently, this test uses a mock EchoData object that only - has attributes. The consolidated flag provided will be used in every to_zarr call (which - is used to write each EchoData group to zarr_path). + Tests to_zarr creation. Currently, this test uses a mock EchoData object that only + has attributes. """ zarr_path = Path("test.zarr") - mock_echodata.to_zarr(str(zarr_path), consolidated=consolidated, overwrite=True) + mock_echodata.to_zarr(str(zarr_path), overwrite=True) + json_path = zarr_path / "zarr.json" - check = True if consolidated else False - zmeta_path = zarr_path / ".zmetadata" - - assert zmeta_path.exists() is check - - if check is True: - check_consolidated(mock_echodata, zmeta_path) + assert json_path.exists() # clean up the zarr file shutil.rmtree(zarr_path) +# TODO: Add test_open_converted with zarr v3 test data since format changed. open_converted works but needs a test. def test_open_converted(ek60_converted_zarr, minio_bucket): # noqa def _check_path(zarr_path): @@ -726,10 +723,10 @@ def test_update_platform_latlon_notimestamp(test_path): ({"time1": 10, "time2": 10}), ], ) -def test_echodata_chunk(chunk_dict): +def test_echodata_chunk(chunk_dict, ek60_path): # Parse Raw File ed = echopype.open_raw( - "echopype/test_data/ek60/DY1801_EK60-D20180211-T164025.raw", sonar_model="EK60" + ek60_path / "DY1801_EK60-D20180211-T164025.raw", sonar_model="EK60" ) # Chunk Echodata object @@ -795,13 +792,13 @@ def test_convert_legacy_versions_ek80(legacy_datatree, legacy_datatree_filename) @pytest.mark.unit -def test_echodata_delete(caplog): +def test_echodata_delete(caplog, ek60_path): """ Check for correct removal behavior and no warnings captured in echodata delete. """ # Open raw using swap file ed = open_raw( - "echopype/test_data/ek60/ncei-wcsd/SH1701/TEST-D20170114-T202932.raw", + ek60_path / "ncei-wcsd/SH1701/TEST-D20170114-T202932.raw", sonar_model="EK60", use_swap=True ) @@ -824,7 +821,7 @@ def test_echodata_delete(caplog): ] if len(zarr_stores) > 0: # Break at the first associated file since there is only one unique file - temp_zarr_path = zarr_stores[0].path + temp_zarr_path = zarr_stores[0].root break if temp_zarr_path: diff --git a/echopype/tests/echodata/utils.py b/echopype/tests/echodata/utils.py index ef2e15aab..e3226e423 100644 --- a/echopype/tests/echodata/utils.py +++ b/echopype/tests/echodata/utils.py @@ -133,18 +133,14 @@ def check_consolidated(echodata: EchoData, zmeta_path: Path) -> None: """ # Check that every group is in # the zmetadata if consolidated - expected_zgroups = [ - os.path.join(p, '.zgroup') if p != 'Top-level' else '.zgroup' - for p in echodata.group_paths - ] + expected_zgroups = echodata.group_paths[1:] with open(zmeta_path) as f: meta_json = json.load(f) file_groups = [ - k - for k in meta_json['metadata'].keys() - if k.endswith('.zgroup') + key for key, value in meta_json['consolidated_metadata']['metadata'].items() + if value.get('node_type') == 'group' ] for g in expected_zgroups: diff --git a/echopype/tests/mask/test_mask.py b/echopype/tests/mask/test_mask.py index 741c31ddb..5b212c78f 100644 --- a/echopype/tests/mask/test_mask.py +++ b/echopype/tests/mask/test_mask.py @@ -16,8 +16,21 @@ _check_freq_diff_source_Sv, ) +# for seafloor +from echopype.mask import detect_seafloor + +# for schoals +from echopype.mask import detect_shoal +from scipy import ndimage as ndi from typing import List, Union, Optional +@pytest.fixture +def ek60_path(test_path): + return test_path["EK60"] + +@pytest.fixture +def ek80_path(test_path): + return test_path["EK80"] def get_mock_freq_diff_data( n: int, @@ -1278,14 +1291,14 @@ def test_apply_mask_channel_variation(source_has_ch, mask, truth_da): ("depth", False), ] ) -def test_apply_mask_dims_using_MVBS(range_var, use_multi_channel_mask): +def test_apply_mask_dims_using_MVBS(range_var, use_multi_channel_mask, ek60_path): """ Check for correct values and dimensions when using `apply_mask` to apply frequency differencing masks to MVBS. """ # Parse Raw File ed = ep.open_raw( - raw_file="echopype/test_data/ek60/DY1801_EK60-D20180211-T164025.raw", + raw_file=ek60_path / "DY1801_EK60-D20180211-T164025.raw", sonar_model="EK60" ) @@ -1358,13 +1371,13 @@ def test_apply_mask_dims_using_MVBS(range_var, use_multi_channel_mask): @pytest.mark.unit -def test_validate_source_ds_and_check_mask_dim_alignment(): +def test_validate_source_ds_and_check_mask_dim_alignment(ek60_path): """ Tests that ValueErrors are raised for `_validate_source_ds_and_check_mask_dim_alignment`. """ # Parse Raw File ed = ep.open_raw( - raw_file="echopype/test_data/ek60/DY1801_EK60-D20180211-T164025.raw", + raw_file=ek60_path / "DY1801_EK60-D20180211-T164025.raw", sonar_model="EK60" ) @@ -1631,4 +1644,365 @@ def test_apply_mask_actual_range_comprehensive(mask_type, test_values, expected_ assert actual_range[0] == actual_min_from_data, \ f"actual_range min doesn't match data min for {test_description}" assert actual_range[1] == actual_max_from_data, \ - f"actual_range max doesn't match data max for {test_description}" \ No newline at end of file + f"actual_range max doesn't match data max for {test_description}" + + +### add test for seafloor + +def _make_ds_Sv_depth(n_ping=100, n_range=40, channel="chan1", dz=1.0, var_name="Sv_corrected"): + """Build a dataset with Sv and a channelized depth grid.""" + ping = np.arange(n_ping) + rs = np.arange(n_range) + depth_1d = rs * dz # depth increases with range_sample + depth_2d = np.tile(depth_1d, (n_ping, 1)) # (ping_time, range_sample) + + ds = xr.Dataset() + ds[var_name] = xr.DataArray( + data=np.full((1, n_ping, n_range), -120.0, dtype=float), + dims=("channel", "ping_time", "range_sample"), + coords={"channel": [channel], "ping_time": ping, "range_sample": rs}, + ) + # make depth have a channel dimension so .sel(channel=...) works + ds["depth"] = xr.DataArray( + data=depth_2d, + dims=("ping_time", "range_sample"), + coords={"ping_time": ping, "range_sample": rs}, + ).expand_dims({"channel": [channel]}).transpose("channel", "ping_time", "range_sample") + + return ds + +@pytest.mark.unit +def test_detect_basic_return(): + ds = _make_ds_Sv_depth() + + # make a sloped bottom band that gets shallower with time + n_ping = ds.sizes["ping_time"] + thickness = 3 + start0 = 22 + slope = -0.02 # bins per ping (negative = shallower with time) + + starts = np.clip( + start0 + np.round(slope * np.arange(n_ping)).astype(int), + 0, + ds.sizes["range_sample"] - thickness + ) + + # fill fake depth with a little noise + rng = np.random.default_rng(42) + sigma = 0.7 + for j, s in enumerate(starts): + noise = rng.normal(0.0, sigma, size=thickness) + ds["Sv_corrected"].isel( + channel=0, ping_time=j, range_sample=slice(s, s + thickness) + )[:] = -35.0 + noise + + # Call dispatcher (basic) + basic_depth = detect_seafloor( + ds, + method="basic", + params={ + "var_name": "Sv_corrected", + "channel": "chan1", + "threshold": (-50.0, -20.0), + "offset_m": 0.3, + "bin_skip_from_surface": 2, + }, + ) + + # Assertions + assert isinstance(basic_depth, xr.DataArray) + assert basic_depth.name == "bottom_depth" + assert basic_depth.dims == ("ping_time",) + assert np.array_equal(basic_depth["ping_time"], ds["ping_time"]) + + idx = xr.DataArray( + starts, + dims=["ping_time"], + coords={"ping_time": ds["ping_time"]} + ) + + # Depth at those first bright bins, one value per ping + depth_at_starts = ds["depth"].isel(range_sample=idx) + + # Expected detector output (subtract the offset in meters) + expected_depth = depth_at_starts.values - 0.3 + + # Compare to the algorithm’s result + assert np.allclose(basic_depth.values, expected_depth, atol=1e-6) + +@pytest.mark.unit +def test_detect_basic_no_detection_defaults_to_bin_skip(): + ds = _make_ds_Sv_depth(n_ping=50, n_range=30, channel="chan1", dz=1.0, var_name="Sv_corrected") + ds["Sv_corrected"][:] = -120.0 # nothing within (-50, -20) => no detection + + out = detect_seafloor( + ds, + method="basic", + params={ + "var_name": "Sv_corrected", + "channel": "chan1", + "threshold": (-50.0, -20.0), + "offset_m": 0.3, + "bin_skip_from_surface": 10, + }, + ) + + assert out.dims == ("ping_time",) + + # Expected = depth at bin 10 (1 sample_range/bin) minus offset, constant across pings + expected_const = ds["depth"].isel(ping_time=0, range_sample=10).item() - 0.3 + expected = np.full(ds.sizes["ping_time"], expected_const) + np.testing.assert_allclose(out.values, expected, atol=1e-6) + + +@pytest.mark.xfail(strict=False, reason="Method may not be implemented yet (#1522)") +def test_detect_seafloor_unknown_method_raises(): + """Expect a ValueError for an unknown method (via dispatcher, if present).""" + from importlib import import_module + + api = import_module("echopype.mask.seafloor_detection.api") # will raise if not present + ds = _make_ds_Sv_depth() + with pytest.raises(ValueError, match="Unsupported.*method"): + api.detect_seafloor( + ds, + method="__bad_method__", + params={"var_name": "Sv_corrected", "channel": "59006-125-2"}, + ) + +@pytest.mark.unit +def test_blackwell_vs_basic_close_local(ek80_path): + """Blackwell vs basic using local test data""" + + raw_path = ek80_path / "ncei-wcsd/SH2306/Hake-D20230811-T165727.raw" + + if not raw_path.is_file(): + pytest.skip(f"Missing EK80 RAW: {raw_path}") + + ed = ep.open_raw(raw_path, sonar_model="EK80") + ds_Sv = ep.calibrate.compute_Sv(ed, waveform_mode="CW", encode_mode="power") + ds_Sv = ep.consolidate.add_depth(ds_Sv, ed, depth_offset=9.8) + + # Pick an available channel + sel_channel = "WBT 400142-15 ES70-7C_ES" + + # attach angles (required by Blackwell) + angle_along = ed["Sonar/Beam_group1"]["angle_alongship"] + angle_athwart = ed["Sonar/Beam_group1"]["angle_athwartship"] + ds_Sv = ds_Sv.assign( + angle_alongship=angle_along, + angle_athwartship=angle_athwart, + ) + + blackwell_depth = detect_seafloor( + ds=ds_Sv, + method="blackwell", + params={ + "channel": sel_channel, + "threshold": [-40, 702, 282], + "offset": 0.3, + "r0": 10, + "r1": 1000, + "wtheta": 28, + "wphi": 52, + }, + ) + + basic_depth = detect_seafloor( + ds_Sv, + method="basic", + params={ + "var_name": "Sv", + "channel": sel_channel, + "threshold": (-40, -20), + "offset_m": 0.3, + "bin_skip_from_surface": 200, + }, + ) + + assert isinstance(blackwell_depth, xr.DataArray) + assert blackwell_depth.dims == ("ping_time",) + assert basic_depth.dims == ("ping_time",) + assert np.array_equal(blackwell_depth["ping_time"], basic_depth["ping_time"]) + + b = blackwell_depth.values + c = basic_depth.values + m = np.isfinite(b) & np.isfinite(c) + assert m.any(), "No overlapping finite detections to compare." + mad = np.median(np.abs(b[m] - c[m])) + assert mad < 5.0, f"Median abs diff too large: {mad:.2f} m" + + +### add test for schoal + +def _make_ds_Sv(n_ping=6, n_range=8, channels=("chan1",)): + coords = {"ping_time": np.arange(n_ping), "range_sample": np.arange(n_range)} + if channels is None: + da = xr.DataArray(np.zeros((n_ping, n_range), float), + dims=("ping_time", "range_sample"), coords=coords, + name="Sv_corrected") + else: + da = xr.DataArray(np.zeros((len(channels), n_ping, n_range), float), + dims=("channel", "ping_time", "range_sample"), + coords={**coords, "channel": list(channels)}, + name="Sv_corrected") + return da.to_dataset() + +@pytest.mark.unit +@pytest.mark.xfail(reason="Unknown-method, raises an error") +def test_detect_shoals_unknown_method_raises(): + ds = _make_ds_Sv(n_ping=2, n_range=2, channels=("59006-125-2",)) + detect_shoal( + ds, + method="__error__", + params={"var_name": "Sv_corrected", "channel": "59006-125-2"}, + ) + +@pytest.mark.unit +def test_weill_basic_gaps_and_sizes(): + """ + Weill: thresholding + vertical/horizontal gap filling in index space. + """ + ds = _make_ds_Sv(n_ping=20, n_range=8, channels=("59006-125-2",)) + + # Background below threshold + ds["Sv_corrected"][:] = -90.0 + + # Vertical pillar at ping 2 with a vertical gap of size 2 (missing 3 and 4) + ds["Sv_corrected"].loc[dict(channel="59006-125-2", ping_time=2, range_sample=[1,2,5,6])] = -60.0 + + # Horizontal bar deeper at range=7 + ds["Sv_corrected"].loc[dict(channel="59006-125-2", ping_time=[0,1,3,4], range_sample=7)] = -55.0 + + # 2x2 at top-left corner, far from both features + ds["Sv_corrected"].loc[dict(channel="59006-125-2", ping_time=[14,15], range_sample=[0,1])] = -57.0 + + # detect + mask = detect_shoal( + ds, + method="weill", + params={ + "var_name": "Sv_corrected", + "channel": "59006-125-2", + "thr": -70, + "maxvgap": 2, + "maxhgap": 1, + "minvlen": 2, + "minhlen": 2, + }, + ) + + # shape/dtype/dims + n_ping = ds.sizes["ping_time"] + n_range = ds.sizes["range_sample"] + + assert mask.dims == ("ping_time", "range_sample") + assert mask.dtype == bool + assert mask.sizes["ping_time"] == n_ping and mask.sizes["range_sample"] == n_range + + # Vertical column at ping=2: expect True from range 1..7 now (row-7 fill bridges at j=2) + col2 = mask.sel(ping_time=2).values + expected_col2 = np.zeros(n_range, dtype=bool) + expected_col2[1:8] = True + np.testing.assert_array_equal(col2, expected_col2) + + # Horizontal row at range=7: expect True at pings 0..4 (with ping 2 filled) + row7 = mask.sel(range_sample=7).values + expected_row7 = np.zeros(n_ping, dtype=bool) + expected_row7[[0, 1, 2, 3, 4]] = True + np.testing.assert_array_equal(row7, expected_row7) + + # The far 2x2 at (ping 14..15, range 0..1) detected with vlen = 2 + blob = mask.sel(ping_time=[14, 15], range_sample=[0, 1]).values + assert blob.any() + + mask = detect_shoal( + ds, + method="weill", + params={ + "var_name": "Sv_corrected", + "channel": "59006-125-2", + "thr": -70, + "maxvgap": 2, + "maxhgap": 1, + "minvlen": 3, # to not detect the 2*2 + "minhlen": 3, # to not detect the 2*2 + }, + ) + + blob = mask.sel(ping_time=[14, 15], range_sample=[0, 1]).values + assert not blob.any() + + +@pytest.mark.unit +def test_weill_dispatcher_and_coords(): + ds = _make_ds_Sv(n_ping=8, n_range=8, channels=("59006-125-2",)) + ds["Sv_corrected"][:] = -90.0 + ds["Sv_corrected"].loc[dict(channel="59006-125-2", ping_time=[1,2], range_sample=[1,2])] = -60.0 + + mask = detect_shoal( + ds, + method="weill", + params={ + "var_name": "Sv_corrected", + "channel": "59006-125-2", + "thr": -70, + "maxvgap": 1, + "maxhgap": 1, + "minvlen": 1, + "minhlen": 1, + }, + ) + + assert mask.dims == ("ping_time", "range_sample") + assert np.array_equal(mask.coords["ping_time"], ds["ping_time"]) + assert np.array_equal(mask.coords["range_sample"], ds["range_sample"]) + + expected = np.zeros((ds.sizes["ping_time"], ds.sizes["range_sample"]), dtype=bool) + expected[1:3, 1:3] = True + np.testing.assert_array_equal(mask.values, expected) + +@pytest.mark.unit +def test_echoview_mincan_no_linking(): + """ + Basic Echoview mask test: create two bright blocks, run detection, check + both are found and remain separate (two connected components). + """ + ds = _make_ds_Sv(n_ping=15, n_range=12, channels=("59006-125-2",)) + ds["Sv_corrected"][:] = -90.0 # background < thr (not kept) + + # Small 2x2 + ds["Sv_corrected"].loc[dict(channel="59006-125-2", ping_time=[1,2], range_sample=[1,2])] = -60.0 + # Large 3x3 moved to range 4..6 (gap at range=3) + ds["Sv_corrected"].loc[dict(channel="59006-125-2", ping_time=[2,3,4], range_sample=[4,5,6])] = -60.0 + + mask = detect_shoal( + ds, + method="echoview", + params={ + "var_name": "Sv_corrected", + "channel": "59006-125-2", + "idim": np.arange(ds.sizes["range_sample"] + 1), # 0..6 (range edges) + "jdim": np.arange(ds.sizes["ping_time"] + 1), # 0..6 (ping edges) + "thr": -70, + "mincan": (2, 2), + "maxlink": (1, 1), + "minsho": (2, 2), + }, + ) + + ####### asserts + assert mask.dims == ("ping_time", "range_sample") + + # 3x3 block must be fully present (pings 2..4, ranges 4..6) + for p in [2, 3, 4]: + for r in [4, 5, 6]: + assert bool(mask.sel(ping_time=p, range_sample=r)), f"Expected True at ({p},{r})" + + # 2x2 block must also be present (pings 1..2, ranges 1..2) + for p in [1, 2]: + for r in [1, 2]: + assert bool(mask.sel(ping_time=p, range_sample=r)), f"Expected True at ({p},{r})" + + # Label the mask and confirm they are separate components (not connected) + _, nlab = ndi.label(mask.values, structure=np.ones((3, 3), dtype=bool)) + assert nlab == 2 diff --git a/echopype/utils/coding.py b/echopype/utils/coding.py index 64e22cdc8..59d1e1ddd 100644 --- a/echopype/utils/coding.py +++ b/echopype/utils/coding.py @@ -3,10 +3,10 @@ import numpy as np import xarray as xr -import zarr from dask.array.core import auto_chunks from dask.utils import parse_bytes from xarray import coding +from zarr.codecs import BloscCodec DEFAULT_TIME_ENCODING = { "units": "nanoseconds since 1970-01-01T00:00:00Z", @@ -16,12 +16,11 @@ COMPRESSION_SETTINGS = { "netcdf4": {"zlib": True, "complevel": 4}, - # zarr compressors were chosen based on xarray results "zarr": { - "float": {"compressor": zarr.Blosc(cname="zstd", clevel=3, shuffle=2)}, - "int": {"compressor": zarr.Blosc(cname="lz4", clevel=5, shuffle=1, blocksize=0)}, - "string": {"compressor": zarr.Blosc(cname="lz4", clevel=5, shuffle=1, blocksize=0)}, - "time": {"compressor": zarr.Blosc(cname="lz4", clevel=5, shuffle=1, blocksize=0)}, + "float": {"compressor": BloscCodec(cname="zstd", clevel=3, shuffle="bitshuffle")}, + "int": {"compressor": BloscCodec(cname="lz4", clevel=5, shuffle="shuffle", blocksize=0)}, + "string": {"compressor": BloscCodec(cname="lz4", clevel=5, shuffle="shuffle", blocksize=0)}, + "time": {"compressor": BloscCodec(cname="lz4", clevel=5, shuffle="shuffle", blocksize=0)}, }, } diff --git a/echopype/utils/io.py b/echopype/utils/io.py index f0f344f87..f311dfe4e 100644 --- a/echopype/utils/io.py +++ b/echopype/utils/io.py @@ -5,6 +5,7 @@ import os import pathlib import platform +import shutil import sys import tempfile import uuid @@ -13,10 +14,10 @@ import fsspec import xarray as xr +import zarr.storage as zs from dask.array import Array as DaskArray from fsspec import AbstractFileSystem, FSMap from fsspec.implementations.local import LocalFileSystem -from zarr.storage import FSStore from ..echodata import EchoData from ..echodata.api import open_converted @@ -69,13 +70,13 @@ def save_file(ds, path, mode, engine, group=None, compression_settings=None, **k # Allows saving both NetCDF and Zarr files from an xarray dataset if engine == "netcdf4": - ds.to_netcdf(path=path, mode=mode, group=group, encoding=encoding, **kwargs) + ds.to_netcdf(path=path, mode=mode, group=group, encoding=encoding, engine=engine, **kwargs) elif engine == "zarr": # Ensure that encoding and chunks match for var, enc in encoding.items(): if isinstance(ds[var].data, DaskArray): ds[var] = ds[var].chunk(enc.get("chunks", {})) - ds.to_zarr(store=path, mode=mode, group=group, encoding=encoding, **kwargs) + ds.to_zarr(store=path.root, mode=mode, group=group, encoding=encoding, **kwargs) else: raise ValueError(f"{engine} is not a supported save format") @@ -473,13 +474,15 @@ def create_temp_zarr_store() -> FSMap: return fsspec.get_mapper(zarr_path) -def delete_zarr_store(store: "FSStore | str", fs: Optional[AbstractFileSystem] = None) -> None: +def delete_zarr_store( + store: "zs.LocalStore | str", fs: Optional[AbstractFileSystem] = None +) -> None: """ Delete the zarr store and all its contents. Parameters ---------- - store : FSStore or str + store : Zarr.Storage or str The store or store path to delete. fs : AbstractFileSystem, optional The fsspec file system to use @@ -498,16 +501,12 @@ def delete_zarr_store(store: "FSStore | str", fs: Optional[AbstractFileSystem] = raise ValueError("Must provide fs if store is a path string") store_path = store else: - # Get the file system, this should already have the - # correct storage options - fs = store.fs - # Get the string path to the store - store_path: str = store.dir_path() + store_path: str = str(store.root) - if fs.exists(store_path): - # Delete the store when it exists - fs.rm(store_path, recursive=True) + if os.path.exists(store_path): + if isinstance(fs, zs.LocalStore): + shutil.rmtree(store_path) # End of utilities for creating temporary swap zarr files ------------------------------ diff --git a/requirements.txt b/requirements.txt index 3fb36a72e..40d4e5533 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,8 +6,8 @@ flox>=0.7.2 fsspec geopy jinja2 -netCDF4>1.6,<1.7.1 -numpy<2 +netCDF4>1.6 +numpy pandas psutil>=5.9.1 pynmea2 @@ -17,4 +17,4 @@ s3fs scipy typing-extensions xarray>=2024.11.0 -zarr>=2,<3 +zarr>=3 diff --git a/setup.cfg b/setup.cfg index 09498df4c..71d5b0329 100644 --- a/setup.cfg +++ b/setup.cfg @@ -15,9 +15,9 @@ classifiers = Operating System :: OS Independent Programming Language :: Python Programming Language :: Python :: 3 :: Only - Programming Language :: Python :: 3.10 Programming Language :: Python :: 3.11 Programming Language :: Python :: 3.12 + Programming Language :: Python :: 3.13 Topic :: Scientific/Engineering author = Wu-Jung Lee author_email = leewujung@gmail.com @@ -29,7 +29,7 @@ platforms = any py_modules = _echopype_version include_package_data = True -python_requires = >=3.9, <3.13 +python_requires = >=3.11, <3.14 setup_requires = setuptools_scm