Exploring Echoregions Regions2D Functionality

Exploring Echoregions Regions2D Functionality#

This notebook parses region values from an Echoview .evr file and creates a region mask for the corresponding Echogram data.

# Importing Packages
import matplotlib.pyplot as plt
import urllib.request
import shutil
import xarray as xr
import numpy as np
from pandas.testing import assert_frame_equal

import echoregions as er

Regions Data Reading#

To start this tutorial, we first download evr data from Echoregions’ Github Repository and parse the .evr file using Echoregions’ read_evr function.

The parsing is based off of the .evr data description shown on Echoview’s website: Region Attributes.

# Set path to notebook data
NOTEBOOK_DATA_PATH = 'https://raw.githubusercontent.com/OSOceanAcoustics/echoregions/store_notebook_data/docs/source/notebook_data'

# Download example EVR File
urllib.request.urlretrieve(f"{NOTEBOOK_DATA_PATH}/transect_notebook.evr","transect_notebook.evr")

# Read EVR file
regions2d = er.read_evr('transect_notebook.evr')

Regions2D as a DataFrame#

regions2d is a specialized object but it has a data attribute which is a simple dataframe.

# Grab regions2d dataframe
regions2d_df = regions2d.data

# Show first 5 rows
regions2d_df
file_name file_type evr_file_format_number echoview_version region_id region_structure_version region_point_count region_selected region_creation_type dummy ... region_bbox_right region_bbox_top region_bbox_bottom region_class region_type region_name time depth region_notes region_detection_settings
0 transect_notebook.evr EVRG 7 13.0.378.44817 1 13 30 0 2 -1 ... 2019-07-02 17:21:19.255 228.020976 295.337409 Hake 1 H20U-Region 0 [2019-07-02T17:16:15.764000000, 2019-07-02T17:... [295.3374093901, 292.1303181193, 290.984928379... [JG, 7/31/19: I removed the shallower part of... []
1 transect_notebook.evr EVRG 7 13.0.378.44817 2 13 32 0 2 -1 ... 2019-07-02 17:16:11.109 68.844398 252.868156 Rockfish 1 Region 1 [2019-07-02T17:15:39.314500000, 2019-07-02T17:... [252.8681563809, 226.5150862912, 203.926740500... [] []
2 transect_notebook.evr EVRG 7 13.0.378.44817 3 13 59 0 2 -1 ... 2019-07-02 17:51:04.562 113.699422 226.917392 Plankton 1 Region 2 [2019-07-02T17:15:07.487500000, 2019-07-02T17:... [113.6994219653, 131.6638187668, 154.113457662... [] []
3 transect_notebook.evr EVRG 7 13.0.378.44817 4 13 11 0 2 -1 ... 2019-07-02 20:26:39.176 232.145737 262.613104 Hake 1 H20U-Region 3 [2019-07-02T20:24:56.339000000, 2019-07-02T20:... [233.7492824627, 259.8641685245, 261.925870055... [JG: AWT20 - continuation of this sign after ... []

4 rows × 22 columns

The regions2d object can be subsetted using the select_region function and with parameters related to region class, time, and depth. For this example let us select a trawl region based on the region_class parameter:

hake_regions = regions2d.select_region(region_class="Hake")

hake_regions
file_name file_type evr_file_format_number echoview_version region_id region_structure_version region_point_count region_selected region_creation_type dummy ... region_bbox_right region_bbox_top region_bbox_bottom region_class region_type region_name time depth region_notes region_detection_settings
0 transect_notebook.evr EVRG 7 13.0.378.44817 1 13 30 0 2 -1 ... 2019-07-02 17:21:19.255 228.020976 295.337409 Hake 1 H20U-Region 0 [2019-07-02T17:16:15.764000000, 2019-07-02T17:... [295.3374093901, 292.1303181193, 290.984928379... [JG, 7/31/19: I removed the shallower part of... []
3 transect_notebook.evr EVRG 7 13.0.378.44817 4 13 11 0 2 -1 ... 2019-07-02 20:26:39.176 232.145737 262.613104 Hake 1 H20U-Region 3 [2019-07-02T20:24:56.339000000, 2019-07-02T20:... [233.7492824627, 259.8641685245, 261.925870055... [JG: AWT20 - continuation of this sign after ... []

2 rows × 22 columns

Now notice that these regions are not closed:

plt.figure(figsize=(20, 6))
for _, point in hake_regions.iterrows():
    plt.plot(point["time"], point["depth"])

plt.xlabel('Ping Time')
plt.ylabel('Depth')
plt.ylim(0, 790)
plt.gca().invert_yaxis()
_images/03dd81ebbb83361ad9b2931cfe9259882eef61cb90e671b9cd569e133a2617d3.png

We can close these regions and re-plot:

hake_regions_closed = regions2d.close_region(region_class="Hake")

plt.figure(figsize=(20, 6))
for _, row in hake_regions_closed.iterrows():
    plt.plot(row["time"], row["depth"])

plt.xlabel('Ping Time')
plt.ylabel('Depth')
plt.ylim(0, 790)
plt.gca().invert_yaxis()
_images/6bc40986676f0ec7f9520545e37bcb5cfbd79d88041fd4d74594f2bf9e6df043.png

To select Hake regions, close regions, and plot regions, we can also just run the following using the object’s plot function:

plt.figure(figsize=(20, 6))
regions2d.plot(region_class="Hake", close_regions=True)
plt.xlabel('Ping Time')
plt.ylabel('Depth')
plt.ylim(0, 790)
plt.gca().invert_yaxis()
_images/6bc40986676f0ec7f9520545e37bcb5cfbd79d88041fd4d74594f2bf9e6df043.png

Echogram Data Reading and Plotting#

Let’s now download a zip file containing backscatter data and plot an echogram created using Echopype. Echopype is a comprehensive software designed for parsing sonar backscatter and conducting scientific computations. It offers a wide array of functionalities, including the plotting of echograms, which are visual representations formed by the echoes received from sonar signals. Additionally, Echopype seamlessly integrates with Echoregions, a specialized software for parsing echogram annotations. For now, we will primarily be working with the ds_Sv["Sv"] backscatter volume data variable, which when plotted becomes an echogram, but there are many more data variables that can be used when working with Echopype.

# Download example Echopype Sv Zarr File
urllib.request.urlretrieve(f"{NOTEBOOK_DATA_PATH}/transect_notebook.zip","transect_notebook.zip")

# Extract the ZIP file
shutil.unpack_archive("transect_notebook.zip", "../../")

# Read volume backscattering strength (Sv) data
ds_Sv = xr.open_dataset('transect_notebook.zarr', engine="zarr")

ds_Sv["Sv"]

<xarray.DataArray 'Sv' (channel: 2, ping_time: 3601, depth: 380)> Size: 22MB
[2736760 values with dtype=float64]
Coordinates:
  * channel    (channel) <U37 296B 'GPT  18 kHz 009072058c8d 1-1 ES18-11' 'GP...
  * ping_time  (ping_time) datetime64[ns] 29kB 2019-07-02T16:00:00 ... 2019-0...
  * depth      (depth) float64 3kB 0.0 2.0 4.0 6.0 ... 752.0 754.0 756.0 758.0
Attributes:
    binning_mode:          physical units
    cell_methods:          ping_time: mean (interval: 5 second comment: ping_...
    long_name:             Mean volume backscattering strength (MVBS, mean Sv...
    ping_time_interval:    5s
    range_meter_interval:  2.0m
    units:                 dB
# Plot the 38 kHz backscatter channel
plt.figure(figsize=(20, 6))
ds_Sv["Sv"].isel(channel=1).plot.pcolormesh(y="depth", yincrease=False, vmin=-70, vmax=-30)

<matplotlib.collections.QuadMesh at 0x772782b90aa0>
_images/a638d2ce16160246020848b156c67e4143acbd11ae8c5e76d79ab054894cb90b.png

Plotting Echogram and Region#

From the two previous plots, we can see how they’re related on both the depth and time dimensions. Now let’s see region annotations overlayed on top of the Echogram dataset.

# Plotting the echogram data and the Hake region
plt.figure(figsize=(20, 6))
for _, point in hake_regions_closed.iterrows():
    plt.plot(point["time"], point["depth"], fillstyle='full', markersize=0.5, color="red")

ds_Sv.Sv.isel(channel=1).T.plot(y="depth", yincrease=False, vmin=-70, vmax=-30)

<matplotlib.collections.QuadMesh at 0x772782c57dd0>
_images/2e7daaad7a57fe211c4f87d550bf78eed959851477c17fb13166f5161a1026c1.png

Creating Echogram Mask with Regions#

From the previous plot, we can see how they’re related on both the depth and time dimensions. Now let’s see a region annotation overlayed on top of the Echogram dataset.

For machine learning training purposes, and more specifically, for computer vision segmentation model training, a region polygon is not enough and so we require masks of the echogram data, which indicate which points are within the region.

The Regions2D class has a function to do this specific task: mask. This mask function is essentially a wrapper over the mask_3D function from regionmask. The 3D in mask_3D is important here since the masks we get out of regions2d will be 3D masks by default (when collapse_to_2d = False in regions2d.mask, which is the default value). This means that the mask will have three dimensions: ping_time, depth, and region_id, where region_id is the unique identifier of a region. The separation of regions via the region_id dimension allows the storage of ‘within region’ and ‘outside region’ information as 1 and 0 respectively.

If the user passes in collapse_to_2d = True they will produce a 2D mask where the mask will have dimensions ping_time, and depth, and the ‘within region’ values will be the region_id values and ‘outside region’ values will be NaN.

Let’s look first at the default collapse_to_2d = False case:

# Use the built in mask function to create a mask
region_mask_ds, region_points = regions2d.region_mask(
    ds_Sv["Sv"].isel(channel=1).drop_vars("channel").compute(),
    region_class="Hake",
    collapse_to_2d= False # region_mask defaults to this so we can comment it out
)

Note that selecting 1 channel here is to ensure that our output is a single channel mask. The number of channels of the data array that is passed into the mask will match the number of channels of the output seafloor mask data array.

The output should be the region mask itself (region_mask_ds) and the points that constitute the region (region_points).

Note that even though we use region_class here to select the region to be masked, what is stored in the dataset as a dimension is the region_id since it is unique and region_class is not:

region_mask_ds

<xarray.Dataset> Size: 22MB
Dimensions:      (region_id: 2, depth: 380, ping_time: 3601)
Coordinates:
  * region_id    (region_id) int64 16B 1 4
  * depth        (depth) float64 3kB 0.0 2.0 4.0 6.0 ... 752.0 754.0 756.0 758.0
  * ping_time    (ping_time) datetime64[ns] 29kB 2019-07-02T16:00:00 ... 2019...
Data variables:
    mask_labels  (region_id) int64 16B 0 1
    mask_3d      (region_id, depth, ping_time) int64 22MB 0 0 0 0 0 ... 0 0 0 0

The values in the region mask should also just be 1s (within the region) and 0s (outside of the region):

print("Unique Values in Region Mask:", np.unique(region_mask_ds["mask_3d"].data))

Unique Values in Region Mask: [0 1]

As a sanity check, let us check that the ping_time and depth dimensions in the region mask region_mask_ds are equal to that in the backscatter data ds_Sv["Sv"]:

print("Region Mask Ping Time Dimension Length:", len(region_mask_ds["ping_time"]))
print("Region Mask Depth Dimension Length:", len(region_mask_ds["depth"]))
print("Echogram Ping Time Dimension Length:", len(ds_Sv["Sv"]["ping_time"]))
print("Echogram Depth Dimension Length:", len(ds_Sv["Sv"]["depth"]))

Region Mask Ping Time Dimension Length: 3601
Region Mask Depth Dimension Length: 380
Echogram Ping Time Dimension Length: 3601
Echogram Depth Dimension Length: 380

A zoomed in plot of the region corresponding to region id 1:

region_mask_ds["mask_3d"].sel(
    region_id=1,
    ping_time=slice("2019-07-02T17:10:00", "2019-07-02T17:25:00")
).plot(y="depth", yincrease=False)

<matplotlib.collections.QuadMesh at 0x77276d3a7e60>
_images/21d39c28ec3d9a69b79c5522f8198868ec02536efa8b3cfe23553fd2bbaf658a.png

A zoomed in plot of the region corresponding to region id 4:

region_mask_ds["mask_3d"].sel(
    region_id=4,
    ping_time=slice("2019-07-02T20:24:00", "2019-07-02T20:27:00")
).plot(y="depth", yincrease=False)

<matplotlib.collections.QuadMesh at 0x77276d07a7e0>
_images/aff17c9290d2d67b0a05f1f09c4e57ac934342091339563f6c42595f7750b393.png

A zoomed in plot of the 38 kHz channel where the mask is 1 and masked region’s region id is 1:

# Get only 38 kHz channel values where the mask is 1 and masked region's region id is 1
mask_exists_Sv = xr.where(
    region_mask_ds["mask_3d"].sel(
        region_id=1,
        drop=False,
    ) == 1,
    ds_Sv["Sv"].isel(channel=1),
    np.nan,
)

# Plot the zoomed in masked Sv
mask_exists_Sv.sel(ping_time=slice("2019-07-02T17:10:00", "2019-07-02T17:25:00")).plot(y="depth", yincrease=False, vmin=-70, vmax=-30)

<matplotlib.collections.QuadMesh at 0x77276d0e12e0>
_images/845208b4520d98f397a1dee0a903268193c2c6d3e1f49e0c2c34d39f8fcebad9.png

Now let us look at what is produced in the collapse_to_2d = True case:

# Set collapse_to_2d = True
region_mask_2d_ds, region_points_2d = regions2d.region_mask(
    ds_Sv["Sv"].isel(channel=1).drop_vars("channel").compute(),
    region_class="Hake",
    collapse_to_2d=True,
)

region_mask_2d_ds["mask_2d"]

<xarray.DataArray 'mask_2d' (depth: 380, ping_time: 3601)> Size: 11MB
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]], shape=(380, 3601))
Coordinates:
  * depth      (depth) float64 3kB 0.0 2.0 4.0 6.0 ... 752.0 754.0 756.0 758.0
  * ping_time  (ping_time) datetime64[ns] 29kB 2019-07-02T16:00:00 ... 2019-0...

The values in the 2d region mask should be 1, 4, and NaN:

print("Unique Values in Region Mask:", np.unique(region_mask_2d_ds["mask_2d"].data))

Unique Values in Region Mask: [ 1.  4. nan]

Note that ds_Sv is the entire dataset containing many data variables that Echoregions does not work with, and so passing in the entire dataset to regions2d.mask will produce an error since it expects a single data array ds_Sv["Sv"]:

# Test ds_Sv
regions2d.region_mask(
    ds_Sv,
    region_class="Hake",
)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[24], line 2
      1 # Test ds_Sv
----> 2 regions2d.region_mask(
      3     ds_Sv,
      4     region_class="Hake",
      5 )

File ~/echoregions/echoregions/regions2d/regions2d.py:597, in Regions2D.region_mask(self, da_Sv, region_id, region_class, mask_name, mask_labels, collapse_to_2d)
    559 """
    560 Mask data from Data Array containing Sv data based off of a Regions2D object
    561 and its regions ids.
   (...)    594     DataFrame containing region_id, depth, and time.
    595 """
    596 if not isinstance(da_Sv, DataArray):
--> 597     raise TypeError(
    598         f"Input da_Sv must be of type DataArray. da_Sv was instead of type {type(da_Sv)}"
    599     )
    601 # Dataframe containing region information.
    602 region_df = self.select_region(region_id, region_class)

TypeError: Input da_Sv must be of type DataArray. da_Sv was instead of type <class 'xarray.core.dataset.Dataset'>

Saving to “.csv” and Reading From “.csv”#

So now that we have our mask and our region polygon points, how do we save them?

We can use the Echoregions read_regions_csv function to first load it onto a region2d object and use the region2d object’s to_csv function to save the regions2d dataframe as a .csv.

# Create new regions2d object
from_mask_regions = er.read_regions_csv(region_points)

# Save to .csv
from_mask_regions.to_csv("from_mask_regions.csv", index=False)

To use er.read_region_csv, the input dataframe must contain (at minimum) columns region_id, depth, and time where each depth entry is a 1-D float array and each time entry is a 1-D datetime64[ns] array:

region_points
region_id depth time
0 1 [295.3374093901, 292.1303181193, 290.984928379... [2019-07-02T17:16:15.764000000, 2019-07-02T17:...
3 4 [233.7492824627, 259.8641685245, 261.925870055... [2019-07-02T20:24:56.339000000, 2019-07-02T20:...

For more information on datetime64[ns] please visit the following: numpy datetime doc.

Now if you need to load this .csv into a regions object we can again use read_regions_csv since it takes in both file locations (Path/str objects) and Pandas DataFrames (as long as the CSV matches the DataFrame format described above):

# Create another new regions2d object
from_csv_regions = er.read_regions_csv("from_mask_regions.csv")

As another sanity check, let’s check if these dataframes are equal (with index reset):

try:
    assert_frame_equal(from_mask_regions.data.reset_index(drop=True), from_csv_regions.data.reset_index(drop=True))
    print("The two DataFrames are equal.")
except AssertionError:
    print("The two DataFrames are not equal.")

The two DataFrames are equal.
Contents