2. EchoMap: Plotting echogram and maps from echo data in the Pacific Hake survey

2.1. Introduction

2.1.1. Goals

• Illustrate how to connect work between echopype and echoshader.

• Illustrate a common workflow for plotting echogram and echomap using echoshader.

2.1.2. Description

This notebook uses EK60 echosounder data collected during the 2017 Joint U.S.-Canada Integrated Ecosystem and Pacific Hake Acoustic Trawl Survey (‘Pacific Hake Survey’) to illustrate a common workflow for data conversion, calibration and visualization using echopype and echoshader.

Ten days of cloud-hosted .raw data files are accessed by echopype directly from an Amazon Web Services (AWS) S3 “bucket” maintained by the NOAA NCEI Water-Column Sonar Data Archive. The total data used are 365 .raw files at approximately 25 MB each (1 Hz pinging rate from first light to dusk), corresponding to approximately 40 GB. With echopype, each file is converted to a standardized representation and saved on Google Drive.

2.1.3. Outline

1. Establish AWS S3 file system connection and Process S3-hosted raw files with echopype.

2. Combine MVBS with coordinate information.

3. Plot echogram, echomap and control widgets with echoshader.

2.1.4. Warning

When Image dimension ping_time is not evenly sampled to relative tolerance of 0.001. There will be an error.

Use echogram.erase_error() to get rid of them.

2.1.5. Note

We encourage importing echopype asep for consistency.

from pathlib import Path
import itertools as it
import datetime as dt
from dateutil import parser as dtparser

import fsspec
from urllib import request

import echopype as ep
import xarray as xr

import panel
import holoviews as hv

from echoshader.echomap import EchoMap

hv.extension("bokeh")

2.2. Establish AWS S3 file system connection and Process S3-hosted raw files with echopype.

base_dpath = Path('./exports')
base_dpath.mkdir(exist_ok=True)

calibrated_dpath = (base_dpath / 'HakeSurvey_10_Days')
calibrated_dpath.mkdir(exist_ok=True)

lon_dpath = (base_dpath / 'Lon_10_Days')
lon_dpath.mkdir(exist_ok=True)

lat_dpath = (base_dpath / 'Lat_10_Days')
lat_dpath.mkdir(exist_ok=True)

fs = fsspec.filesystem('s3', anon=True)

bucket = "ncei-wcsd-archive"
rawdirpath = "data/raw/Bell_M._Shimada/SH1707/EK60"

s3rawfiles = fs.glob(f"{bucket}/{rawdirpath}/*.raw")

start_datetime = dt.datetime(2017, 7, 24, 0, 0)
end_datetime = dt.datetime(2017, 8, 2, 23, 59)

date_list = []

for x in range(0, (end_datetime-start_datetime).days + 1):
    date = start_datetime + dt.timedelta(days = x)
    date_list.append(str(date.month).zfill(2)+str(date.day).zfill(2))

s3rawfiles = [
    s3path for s3path in s3rawfiles 
    if any([f"D2017{datestr}" in s3path for datestr in date_list])
]

print(f"There are {len(s3rawfiles)} target raw files available")

for s3rawfpath in s3rawfiles:
    raw_fpath = Path(s3rawfpath)
    try:
        # Access file directly from S3 to create a converted EchoData object in memory
        ed = ep.open_raw(
            f"s3://{s3rawfpath}",
            sonar_model='EK60',
            storage_options={'anon': True}
        )
        
        
        # Use the EchoData object "ed" to generate calibrated and
        # computed MVBS files that will be saved to netcdf
        ds_Sv = ep.calibrate.compute_Sv(ed)
        
        ds_MVBS = ep.preprocess.compute_MVBS(
            ds_Sv,
            range_meter_bin=5,  # in meters
            ping_time_bin='1800s'  # in seconds
        )
        
        ds_MVBS.to_netcdf(calibrated_dpath / f"MVBS_{raw_fpath.stem}.nc")

        ds_lon = ed['Platform'].longitude
        
        ds_lon.to_netcdf(lon_dpath / f"MVBS_{raw_fpath.stem}.nc")
        
        ds_lat = ed['Platform'].latitude
        
        ds_lat.to_netcdf(lat_dpath / f"MVBS_{raw_fpath.stem}.nc")
        
    except Exception as e:
        print(f"Failed to process raw file {raw_fpath.name}: {e}")

2.3 Combine MVBS with coordinate information.

MVBS_ds = xr.open_mfdataset(
    str(calibrated_dpath / 'MVBS_*.nc'), 
    data_vars='minimal', coords='minimal',
    combine='by_coords'
)

longitude = xr.open_mfdataset(
    str(lon_dpath / '*.nc'),
    data_vars='minimal', coords='minimal',
    combine='by_coords'
)

latitude = xr.open_mfdataset(
    str(lat_dpath / '*.nc'),
    data_vars='minimal', coords='minimal',
    combine='by_coords'
)

lon = longitude["longitude"]

lat = latitude["latitude"]

lon=lon.interp(time1=MVBS_ds["ping_time"])
lat=lat.interp(time1=MVBS_ds["ping_time"])

MVBS_ds["longitude"]=lon
MVBS_ds["latitude"]=lat

import datetime
history = (
        f"{datetime.datetime.utcnow()} +00:00. "
        "Interpolated from Platform latitude/longitude."
    )
MVBS_ds["latitude"] = MVBS_ds["latitude"].assign_attrs({"history": history})
MVBS_ds["longitude"] = MVBS_ds["longitude"].assign_attrs({"history": history})

MVBS_ds

2.4 Plot echogram, echomap and control widgets with echoshader.

# Calibratd data is stored in Google Drive
url = 'https://drive.google.com/uc?export=download&id=1E5mdgALsPApD2-vFoetIZXKyo_y3NUb9'

def urllib_download():
    request.urlretrieve(url, 'test_EchoMap.nc')

urllib_download()   

MVBS_ds = xr.open_mfdataset(
    paths = 'test_EchoMap.nc', 
    data_vars = 'minimal', coords='minimal',
    combine = 'by_coords'
)

MVBS_ds

<xarray.Dataset>
Dimensions:            (channel: 4, ping_time: 875, echo_range: 150)
Coordinates:
  * channel            (channel) object 'GPT  18 kHz 009072058c8d 1-1 ES18-11...
  * ping_time          (ping_time) datetime64[ns] 2017-07-24T19:30:00 ... 201...
  * echo_range         (echo_range) float64 0.0 5.0 10.0 ... 735.0 740.0 745.0
    time1              (ping_time) datetime64[ns] dask.array<chunksize=(875,), meta=np.ndarray>
Data variables:
    Sv                 (channel, ping_time, echo_range) float64 dask.array<chunksize=(4, 875, 150), meta=np.ndarray>
    frequency_nominal  (channel) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    longitude          (ping_time) float64 dask.array<chunksize=(875,), meta=np.ndarray>
    latitude           (ping_time) float64 dask.array<chunksize=(875,), meta=np.ndarray>
Attributes:
    processing_software_name:     echopype
    processing_software_version:  0.6.0
    processing_time:              2022-08-20T07:11:30Z
    processing_function:          preprocess.compute_MVBS

xarray.Dataset

• Dimensions:
  • channel: 4
  • ping_time: 875
  • echo_range: 150
• Coordinates: (4)
  • channel
    (channel)
    object
    'GPT 18 kHz 009072058c8d 1-1 ES...

    long_name :

    Vendor channel ID

    array(['GPT  18 kHz 009072058c8d 1-1 ES18-11',
           'GPT  38 kHz 009072058146 2-1 ES38B',
           'GPT 120 kHz 00907205a6d0 4-1 ES120-7C',
           'GPT 200 kHz 00907205aec4 3-1 ES200-7C'], dtype=object)

  • ping_time
    (ping_time)
    datetime64[ns]
    2017-07-24T19:30:00 ... 2017-08-03

    long_name :

    Ping time

    standard_name :

    time

    axis :

    T

    array(['2017-07-24T19:30:00.000000000', '2017-07-24T20:00:00.000000000',
           '2017-07-24T20:00:00.000000000', ..., '2017-08-02T23:30:00.000000000',
           '2017-08-02T23:30:00.000000000', '2017-08-03T00:00:00.000000000'],
          dtype='datetime64[ns]')

  • echo_range
    (echo_range)
    float64
    0.0 5.0 10.0 ... 735.0 740.0 745.0

    long_name :

    Range distance

    units :

    m

    array([  0.,   5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.,  50.,  55.,
            60.,  65.,  70.,  75.,  80.,  85.,  90.,  95., 100., 105., 110., 115.,
           120., 125., 130., 135., 140., 145., 150., 155., 160., 165., 170., 175.,
           180., 185., 190., 195., 200., 205., 210., 215., 220., 225., 230., 235.,
           240., 245., 250., 255., 260., 265., 270., 275., 280., 285., 290., 295.,
           300., 305., 310., 315., 320., 325., 330., 335., 340., 345., 350., 355.,
           360., 365., 370., 375., 380., 385., 390., 395., 400., 405., 410., 415.,
           420., 425., 430., 435., 440., 445., 450., 455., 460., 465., 470., 475.,
           480., 485., 490., 495., 500., 505., 510., 515., 520., 525., 530., 535.,
           540., 545., 550., 555., 560., 565., 570., 575., 580., 585., 590., 595.,
           600., 605., 610., 615., 620., 625., 630., 635., 640., 645., 650., 655.,
           660., 665., 670., 675., 680., 685., 690., 695., 700., 705., 710., 715.,
           720., 725., 730., 735., 740., 745.])

  • time1
    (ping_time)
    datetime64[ns]
    dask.array<chunksize=(875,), meta=np.ndarray>

    long_name :

    Ping time

    standard_name :

    time

    axis :

    T

    Array Chunk Bytes 6.84 kiB 6.84 kiB Shape (875,) (875,) Count 2 Graph Layers 1 Chunks Type datetime64[ns] numpy.ndarray

• Data variables: (4)
  • Sv
    (channel, ping_time, echo_range)
    float64
    dask.array<chunksize=(4, 875, 150), meta=np.ndarray>

    long_name :

    Mean volume backscattering strength (MVBS, mean Sv re 1 m-1)

    units :

    dB

    actual_range :

    [-95.33 14.31]

    cell_methods :

    ping_time: mean (interval: 30 minute comment: ping_time is the interval start) echo_range: mean (interval: 5 meter comment: echo_range is the interval start)

    binning_mode :

    physical units

    range_meter_interval :

    5m

    ping_time_interval :

    1800s

    Array Chunk Bytes 4.01 MiB 4.01 MiB Shape (4, 875, 150) (4, 875, 150) Count 2 Graph Layers 1 Chunks Type float64 numpy.ndarray

  • frequency_nominal
    (channel)
    float64
    dask.array<chunksize=(4,), meta=np.ndarray>

    units :

    Hz

    long_name :

    Transducer frequency

    valid_min :

    0.0

    standard_name :

    sound_frequency

    Array Chunk Bytes 32 B 32 B Shape (4,) (4,) Count 2 Graph Layers 1 Chunks Type float64 numpy.ndarray

  • longitude
    (ping_time)
    float64
    dask.array<chunksize=(875,), meta=np.ndarray>

    long_name :

    Platform longitude

    standard_name :

    longitude

    units :

    degrees_east

    valid_range :

    (-180.0, 180.0)

    history :

    2022-08-27 02:58:34.500227 +00:00. Interpolated from Platform latitude/longitude.

    Array Chunk Bytes 6.84 kiB 6.84 kiB Shape (875,) (875,) Count 2 Graph Layers 1 Chunks Type float64 numpy.ndarray

  • latitude
    (ping_time)
    float64
    dask.array<chunksize=(875,), meta=np.ndarray>

    long_name :

    Platform latitude

    standard_name :

    latitude

    units :

    degrees_north

    valid_range :

    (-90.0, 90.0)

    history :

    2022-08-27 02:58:34.500227 +00:00. Interpolated from Platform latitude/longitude.

    Array Chunk Bytes 6.84 kiB 6.84 kiB Shape (875,) (875,) Count 2 Graph Layers 1 Chunks Type float64 numpy.ndarray

• Attributes: (4)

  processing_software_name :

  echopype

  processing_software_version :

  0.6.0

  processing_time :

  2022-08-20T07:11:30Z

  processing_function :

  preprocess.compute_MVBS

There using EchoMap, wchich is a subclass of Echogram.

echomap = EchoMap(MVBS_ds)

Below shows how to get map in select box.

%%capture --no-display
panel.Row(echomap.widgets , panel.Column(echomap.view_gram, echomap.view_box_map))

WARNING:param.Image01534: Image dimension ping_time is  not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.
WARNING:param.Image01534: Image dimension ping_time is  not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.

Below shows how to get map in select box and map responding to full-screen echogram.

panel.Row(panel.Column(echomap.view_all_map))

Below shows how to get curtain in select box.

panel.Row(echomap.view_curtain)

Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process.