Bathymetry#

Elevation relative to sea level / Sea-floor height (above Lowest Astronomical Tide datum)

import os
from zipfile import BadZipFile, ZipFile

import cartopy.crs as ccrs
import contextily as cx
import matplotlib.pyplot as plt
import numpy as np
import rasterio as rio
import seaborn as sns
import xarray as xr
from matplotlib_scalebar.scalebar import ScaleBar

from h2ss import data as rd

plt.rcParams["xtick.major.size"] = 0
plt.rcParams["ytick.major.size"] = 0

# base data download directory
DATA_DIR = os.path.join("data", "bathymetry")

# DTM tile D4
FILE_NAME = "D4_2022.nc.zip"

URL = f"https://downloads.emodnet-bathymetry.eu/v11/{FILE_NAME}"

DATA_FILE = os.path.join(DATA_DIR, FILE_NAME)

# basemap cache directory
cx.set_cache_dir(os.path.join("data", "basemaps"))

rd.download_data(url=URL, data_dir=DATA_DIR, file_name=FILE_NAME)

Data 'D4_2022.nc.zip' already exists in 'data/bathymetry'.
Data downloaded on: 2024-01-06 16:28:46.418620+00:00
Download URL: https://downloads.emodnet-bathymetry.eu/v11/D4_2022.nc.zip
SHA256 hash: 8e6fb771887804dc096e18887ec434505a08ea5092981a33065afdf3225780e0

ZipFile(DATA_FILE).namelist()

['D4_2022.nc']

# extract the archive
try:
    z = ZipFile(DATA_FILE)
    z.extractall(DATA_DIR)
except BadZipFile:
    print("There were issues with the file", DATA_FILE)

data = xr.open_dataset(
    os.path.join(DATA_DIR, "D4_2022.nc"), decode_coords="all"
)

data = data.chunk({"lat": 1000, "lon": 1000, "cdi_index_count": 1000})

data

<xarray.Dataset> Size: 3GB
Dimensions:             (lon: 9484, lat: 9004, cdi_index_count: 21797)
Coordinates:
  * lon                 (lon) float64 76kB -6.377 -6.376 -6.374 ... 3.501 3.502
  * lat                 (lat) float64 72kB 52.5 52.5 52.5 ... 61.87 61.88 61.88
    crs                 float64 8B ...
Dimensions without coordinates: cdi_index_count
Data variables:
    elevation           (lat, lon) float32 342MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    value_count         (lat, lon) float64 683MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    cdi_index           (lat, lon) float64 683MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    interpolation_flag  (lat, lon) float32 342MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    elevation_max       (lat, lon) float32 342MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    elevation_min       (lat, lon) float32 342MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    stdev               (lat, lon) float32 342MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    cdi_reference       (cdi_index_count) object 174kB dask.array<chunksize=(1000,), meta=np.ndarray>
Attributes:
    dtm_convention_version:  1.0
    Conventions:             SeaDataNet_1.0 CF1.6
    title:                   The EMODnet Grid
    institution:             On behalf of the EMODnet project, http://www.emo...
    source:                  source of the data can be found in the dataset o...
    references:              WORK IN PROGRESS 2020 lastest release is DOI: 10...
    comment:                 The data in the EMODnet Grid should not be used ...
    history:                 NetCDF file created with GGSgc NC_Makegrid versi...
data.rio.crs

CRS.from_wkt('GEOGCS["undefined",DATUM["undefined",SPHEROID["undefined",6378137,298.257223563]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]')

data.rio.bounds()

(-6.377083333333333, 52.49791666666667, 3.502083333332933, 61.877083333324805)

data.rio.resolution()

(0.0010416666666666244, 0.0010416666666657193)

# read Kish Basin data and extent
ds, extent = rd.read_dat_file(dat_path=os.path.join("data", "kish-basin"))

xmin, ymin, xmax, ymax = extent.total_bounds

shape = rd.halite_shape(dat_xr=ds).buffer(1000).buffer(-1000)

# reproject bathymetry data to the Kish Basin data's CRS
data_ = data.rio.reproject(rd.CRS).rio.clip(extent)

data_.rio.resolution()

(92.06798784599279, -92.06798784599353)

def plot_bath_map(xds, levels=15, cmap="mako", vmax=None, vmin=None):
    """Plotting helper function"""
    plt.figure(figsize=(10, 7))
    ax = plt.axes(projection=ccrs.epsg(rd.CRS))
    xds.plot.contourf(
        ax=ax,
        robust=True,
        cmap=cmap,
        levels=levels,
        vmax=vmax,
        vmin=vmin,
        extend="both",
    )
    shape.boundary.plot(ax=ax, color="white")

    plt.ylim(ymin - 3000, ymax + 3000)
    plt.xlim(xmin - 3000, xmax + 3000)

    cx.add_basemap(
        ax, source=cx.providers.CartoDB.Positron, crs=rd.CRS, zoom=10
    )
    ax.gridlines(
        draw_labels={"bottom": "x", "left": "y"},
        alpha=0.25,
        color="darkslategrey",
    )
    ax.add_artist(
        ScaleBar(1, box_alpha=0, location="lower right", color="darkslategrey")
    )

    plt.title(None)
    plt.tight_layout()
    plt.show()

plot_bath_map(data_["elevation"])
../_images/20553e15152d947d55ac6cff4f4c9fefa301124103db111e4d2bf4a493ac8c8d.png

Reproject bathymetry to match the resolution of the Kish Basin data#

data_ = data.rename({"lon": "x", "lat": "y"}).rio.reproject_match(
    ds, resampling=rio.enums.Resampling.bilinear
)

data_

<xarray.Dataset> Size: 2MB
Dimensions:             (x: 218, y: 237, cdi_index_count: 21797)
Coordinates:
    crs                 int64 8B 0
  * x                   (x) float64 2kB 6.966e+05 6.968e+05 ... 7.4e+05
  * y                   (y) float64 2kB 5.936e+06 5.936e+06 ... 5.889e+06
Dimensions without coordinates: cdi_index_count
Data variables:
    elevation           (y, x) float32 207kB -15.55 -16.94 ... -60.4 -60.26
    value_count         (y, x) float64 413kB 945.9 1.758e+03 ... 1.0 1.0
    cdi_index           (y, x) float64 413kB 1.856e+04 1.856e+04 ... 1.473e+04
    interpolation_flag  (y, x) float32 207kB 0.0 0.0 0.0 0.0 ... nan nan nan nan
    elevation_max       (y, x) float32 207kB -16.18 -17.76 ... -60.4 -60.26
    elevation_min       (y, x) float32 207kB -15.08 -16.36 ... -60.4 -60.26
    stdev               (y, x) float32 207kB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    cdi_reference       (cdi_index_count) object 174kB dask.array<chunksize=(1000,), meta=np.ndarray>
Attributes:
    dtm_convention_version:  1.0
    Conventions:             SeaDataNet_1.0 CF1.6
    title:                   The EMODnet Grid
    institution:             On behalf of the EMODnet project, http://www.emo...
    source:                  source of the data can be found in the dataset o...
    references:              WORK IN PROGRESS 2020 lastest release is DOI: 10...
    comment:                 The data in the EMODnet Grid should not be used ...
    history:                 NetCDF file created with GGSgc NC_Makegrid versi...
data_.rio.crs

CRS.from_epsg(23029)

data_.rio.bounds()

(696500.0, 5889100.0, 740100.0, 5936500.0)

data_.rio.resolution()

(200.0, -200.0)

val = data_["elevation"].values.flatten()
min(val[~np.isnan(val)])

-166.57455

max(val[~np.isnan(val)])

-0.90690255

plot_bath_map(data_["elevation"], levels=14, vmax=0, vmin=-130)
../_images/60d9e2d87ce580edcb2f57c53a30d32f06d5d04b022a432fd0a4ec854f11066b.png

Adjust Kish Basin depth from sea level to seabed#

ds = ds.assign(TopDepthSeabed=ds["TopDepth"] + data_["elevation"])
ds = ds.assign(BaseDepthSeabed=ds["BaseDepth"] + data_["elevation"])

# ds["TopDepthSeabed"].attrs = ds["TopDepth"].attrs
# ds["BaseDepthSeabed"].attrs = ds["BaseDepth"].attrs

ds

<xarray.Dataset> Size: 10MB
Dimensions:          (halite: 4, y: 237, x: 218)
Coordinates:
  * y                (y) float64 2kB 5.936e+06 5.936e+06 ... 5.889e+06 5.889e+06
  * x                (x) float64 2kB 6.966e+05 6.968e+05 ... 7.398e+05 7.4e+05
    spatial_ref      int64 8B 0
  * halite           (halite) <U8 128B 'Fylde' 'Mythop' 'Preesall' 'Rossall'
    crs              int64 8B 0
Data variables:
    BaseDepth        (halite, y, x) float64 2MB nan nan nan nan ... nan nan nan
    Thickness        (halite, y, x) float64 2MB nan nan nan nan ... nan nan nan
    TopDepth         (halite, y, x) float64 2MB nan nan nan nan ... nan nan nan
    TopTWT           (halite, y, x) float64 2MB nan nan nan nan ... nan nan nan
    TopDepthSeabed   (halite, y, x) float64 2MB nan nan nan nan ... nan nan nan
    BaseDepthSeabed  (halite, y, x) float64 2MB nan nan nan nan ... nan nan nan
val = ds["TopDepthSeabed"].values.flatten()
min(val[~np.isnan(val)])

69.23700712890626

max(val[~np.isnan(val)])

3234.9487763793945

val = ds["BaseDepthSeabed"].values.flatten()
min(val[~np.isnan(val)])

73.09180877075195

max(val[~np.isnan(val)])

3475.1658741638184

plot_bath_map(
    ds["TopDepthSeabed"].sel(halite="Rossall"), cmap="jet", vmin=150, vmax=2250
)
../_images/1fd806684704ee673d9a155c9c52a03e2caac482540f9a00cb6d8eaf024a6cc8.png 
plot_bath_map(ds["BaseDepthSeabed"].sel(halite="Rossall"), cmap="jet")
../_images/9aaa1ba32df5ee58c8d6114d04ac0756ad26b806c9abf612d05837c00f80ad2c.png 
bath = rd.bathymetry_layer(
    dat_extent=extent.buffer(3000),
    bathymetry_path=os.path.join("data", "bathymetry"),
)

plt.figure(figsize=(10, 7))
ax2 = plt.axes(projection=ccrs.epsg(rd.CRS))
bath["elevation"].plot.contourf(cmap="mako", levels=10, robust=True)
CS = bath["elevation"].plot.contour(
    # colors="darkslategrey",
    colors="white",
    levels=10,
    linewidths=0.5,
    linestyles="solid",
    # alpha=0.5,
    robust=True,
)
plt.clabel(CS, inline=True)
cx.add_basemap(ax2, source=cx.providers.CartoDB.VoyagerNoLabels, crs=rd.CRS)
cx.add_basemap(ax2, source=cx.providers.CartoDB.VoyagerOnlyLabels, crs=rd.CRS)
ax2.gridlines(
    draw_labels={"bottom": "x", "left": "y"},
    alpha=0.25,
    color="darkslategrey",
)
plt.title(None)
plt.tight_layout()
plt.show()
../_images/116127fbddfad6537aee4e9f3424934c0d7496a8bd91f9f990e0818fe28736f7.png
On this page
Show Source