Satellite Derived Bathymetry - Caballero & Stumpf Method¶

This example uses Caballero & Stumpf method to derive bathymetry as mentioned in this article: https://doi.org/10.3390/rs12030451

In [1]:

import matplotlib.pyplot as plt
import contextily as ctx
import numpy as np
import cmocean

from sensingpy import reader, plot
from sensingpy.bathymetry.models import stumpf_pseudomodel, LinearModel
from sensingpy.bathymetry.metrics import ValidationSummary
from sensingpy.bathymetry.plot import CalibrationPlot, ValidationPlot, get_approximate_density
from sensingpy.masks import is_valid
from sensingpy.preprocessing.outliers import IQR
from sensingpy.enums import SENTINEL2_BANDS
from sensingpy.selector import composite_indices
from sensingpy.image import compose

import matplotlib.pyplot as plt import contextily as ctx import numpy as np import cmocean from sensingpy import reader, plot from sensingpy.bathymetry.models import stumpf_pseudomodel, LinearModel from sensingpy.bathymetry.metrics import ValidationSummary from sensingpy.bathymetry.plot import CalibrationPlot, ValidationPlot, get_approximate_density from sensingpy.masks import is_valid from sensingpy.preprocessing.outliers import IQR from sensingpy.enums import SENTINEL2_BANDS from sensingpy.selector import composite_indices from sensingpy.image import compose

In [2]:

import warnings

warnings.filterwarnings("ignore")

import warnings warnings.filterwarnings("ignore")

Load Satellite Images¶

In [3]:

P_SDB_GREEN = 'pSDB Green'
DEPTH = 'DEPTH'

P_SDB_GREEN = 'pSDB Green' DEPTH = 'DEPTH'

In [4]:

filenames = [
    r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181022_S2B_AGV_EXT_29SPA_L2W.nc",
    r"E:\batimetria\algarve\satelite\local_original\Pleamar\20180919_S2B_AGV_EXT_29SPA_L2W.nc",
    r"E:\batimetria\algarve\satelite\local_original\Pleamar\20180922_S2B_AGV_EXT_29SPA_L2W.nc",
    r"E:\batimetria\algarve\satelite\local_original\Pleamar\20180924_S2A_AGV_EXT_29SPA_L2W.nc",
    r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181004_S2A_AGV_EXT_29SPA_L2W.nc",
    r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181007_S2A_AGV_EXT_29SPA_L2W.nc",
    r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181017_S2A_AGV_EXT_29SPA_L2W.nc",
    r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181019_S2B_AGV_EXT_29SPA_L2W.nc",
]

images = []
for filename in filenames:
    image = reader.open(filename)
    image.rename_by_enum(SENTINEL2_BANDS)
    image.rename({'Rrs_B2' : 'Blue', 'Rrs_B3' : 'Green', 'Rrs_B4' : 'Red'})
    image.keep_bands( ['Blue', 'Green', 'Red'] )

    # pSDB Green
    image[P_SDB_GREEN] = IQR( stumpf_pseudomodel(image['Blue'], image['Green']), distance = 1.5 )
    
    images.append(image)

filenames = [ r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181022_S2B_AGV_EXT_29SPA_L2W.nc", r"E:\batimetria\algarve\satelite\local_original\Pleamar\20180919_S2B_AGV_EXT_29SPA_L2W.nc", r"E:\batimetria\algarve\satelite\local_original\Pleamar\20180922_S2B_AGV_EXT_29SPA_L2W.nc", r"E:\batimetria\algarve\satelite\local_original\Pleamar\20180924_S2A_AGV_EXT_29SPA_L2W.nc", r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181004_S2A_AGV_EXT_29SPA_L2W.nc", r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181007_S2A_AGV_EXT_29SPA_L2W.nc", r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181017_S2A_AGV_EXT_29SPA_L2W.nc", r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181019_S2B_AGV_EXT_29SPA_L2W.nc", ] images = [] for filename in filenames: image = reader.open(filename) image.rename_by_enum(SENTINEL2_BANDS) image.rename({'Rrs_B2' : 'Blue', 'Rrs_B3' : 'Green', 'Rrs_B4' : 'Red'}) image.keep_bands( ['Blue', 'Green', 'Red'] ) # pSDB Green image[P_SDB_GREEN] = IQR( stumpf_pseudomodel(image['Blue'], image['Green']), distance = 1.5 ) images.append(image)

Compute clearest image¶

In [5]:

indices = composite_indices( np.array( [ image[P_SDB_GREEN] for image in images ] ), method = np.nanargmax )
multi_image = compose(images, indices)

indices = composite_indices( np.array( [ image[P_SDB_GREEN] for image in images ] ), method = np.nanargmax ) multi_image = compose(images, indices)

Load In-Situ Depths¶

In [6]:

in_situ = reader.open(r"E:\batimetria\algarve\in_situ\scaled\Armona\2018_10m.tif")
in_situ.dropna()

fig, ax = plot.get_geofigure(in_situ.crs, 1, 1)
ax, mappable = plot.plot_band(in_situ, DEPTH, ax = ax, cmap = cmocean.cm.deep)
ctx.add_basemap(ax, crs = in_situ.crs, attribution=False)
fig.colorbar(mappable, label = 'Depth (m)')
plot.add_gridlines(ax)

in_situ = reader.open(r"E:\batimetria\algarve\in_situ\scaled\Armona\2018_10m.tif") in_situ.dropna() fig, ax = plot.get_geofigure(in_situ.crs, 1, 1) ax, mappable = plot.plot_band(in_situ, DEPTH, ax = ax, cmap = cmocean.cm.deep) ctx.add_basemap(ax, crs = in_situ.crs, attribution=False) fig.colorbar(mappable, label = 'Depth (m)') plot.add_gridlines(ax)

Out[6]:

(<GeoAxes: >, <cartopy.mpl.gridliner.Gridliner at 0x22e81f87590>)

Align In-Situ with Satellite dimensions¶

In [7]:

in_situ.align(multi_image)
in_situ.dropna()
multi_image_clip = multi_image.clip([in_situ.bbox], inplace = False)

in_situ.align(multi_image) in_situ.dropna() multi_image_clip = multi_image.clip([in_situ.bbox], inplace = False)

Impact of multi-image over pSDB green¶

In [8]:

for image in images + [multi_image]:
    data = image.clip([in_situ.bbox], inplace = False)[P_SDB_GREEN]

    x, y = data.ravel(), in_situ[DEPTH].ravel()
    valid = is_valid(x) & is_valid(y)
    x, y = x[valid], y[valid]

    fig, ax = plt.subplots(1, 1)
    x, y, z, norm = get_approximate_density(x, y, bins = 100)
    ax.scatter(x, y, c = z, norm=norm, cmap = 'viridis_r', s = 0.1)
    ax.set_xlabel(f'{P_SDB_GREEN} (DII)')
    ax.set_ylabel(f'depth (m)')
    ax.set_xlim(0.5, 1.5)

for image in images + [multi_image]: data = image.clip([in_situ.bbox], inplace = False)[P_SDB_GREEN] x, y = data.ravel(), in_situ[DEPTH].ravel() valid = is_valid(x) & is_valid(y) x, y = x[valid], y[valid] fig, ax = plt.subplots(1, 1) x, y, z, norm = get_approximate_density(x, y, bins = 100) ax.scatter(x, y, c = z, norm=norm, cmap = 'viridis_r', s = 0.1) ax.set_xlabel(f'{P_SDB_GREEN} (DII)') ax.set_ylabel(f'depth (m)') ax.set_xlim(0.5, 1.5)

In-situ and Pseudomodel¶

In [9]:

fig, axs = plot.get_geofigure(multi_image_clip.crs, 1, 2)
plot.plot_band(multi_image_clip, P_SDB_GREEN, axs[0], cmap = cmocean.cm.deep)
_, mappable = plot.plot_band(in_situ, DEPTH, ax = axs[1], cmap = cmocean.cm.deep)

ctx.add_basemap(axs[0], crs = multi_image_clip.crs, attribution=False)
ctx.add_basemap(axs[1], crs = in_situ.crs, attribution=False)


axs[0].set_title('pSDB Green')
axs[1].set_title('In-Situ Depth')
plot.add_gridlines(axs[0])
plot.add_gridlines(axs[1])

fig, axs = plot.get_geofigure(multi_image_clip.crs, 1, 2) plot.plot_band(multi_image_clip, P_SDB_GREEN, axs[0], cmap = cmocean.cm.deep) _, mappable = plot.plot_band(in_situ, DEPTH, ax = axs[1], cmap = cmocean.cm.deep) ctx.add_basemap(axs[0], crs = multi_image_clip.crs, attribution=False) ctx.add_basemap(axs[1], crs = in_situ.crs, attribution=False) axs[0].set_title('pSDB Green') axs[1].set_title('In-Situ Depth') plot.add_gridlines(axs[0]) plot.add_gridlines(axs[1])

Out[9]:

(<GeoAxes: title={'center': 'In-Situ Depth'}>,
 <cartopy.mpl.gridliner.Gridliner at 0x22e94917150>)

Calibration¶

Select random depths¶

arginterval_choice returns a list of indexes. the band is transformed to 1D internaly, so we must transform our bands to 1D when masking.

In [10]:

selection = in_situ.arginterval_choice(band = DEPTH, size = 3, intervals = range(0, 10))

selection = in_situ.arginterval_choice(band = DEPTH, size = 3, intervals = range(0, 10))

Extract depths and pseudomodel as 1D-arrays¶

In [11]:

depths = in_situ.select(DEPTH).ravel()
p_green = multi_image_clip.select(P_SDB_GREEN).ravel()

depths = in_situ.select(DEPTH).ravel() p_green = multi_image_clip.select(P_SDB_GREEN).ravel()

Fit model¶

In [12]:

to_cal_depths = depths[selection]
to_cal_p_greens = p_green[selection]
no_nans = is_valid(to_cal_depths) & is_valid(to_cal_p_greens)

calibration = LinearModel().fit(to_cal_p_greens[no_nans], to_cal_depths[no_nans])

print(f'N: {len(to_cal_depths[no_nans])}, Calibration: {calibration}')

to_cal_depths = depths[selection] to_cal_p_greens = p_green[selection] no_nans = is_valid(to_cal_depths) & is_valid(to_cal_p_greens) calibration = LinearModel().fit(to_cal_p_greens[no_nans], to_cal_depths[no_nans]) print(f'N: {len(to_cal_depths[no_nans])}, Calibration: {calibration}')

N: 27, Calibration: R: 0.8965 | y = 59.550x-55.668

Plot model¶

In [13]:

fig, ax = plt.subplots(1, 1, figsize = (6, 6))
cal_plot = CalibrationPlot(legend_font_size = 15)
cal_plot.add_calibration_scatter(calibration, to_cal_p_greens[no_nans], to_cal_depths[no_nans], ax = ax)
cal_plot.add_labels(ax, title = f'{P_SDB_GREEN} vs {DEPTH}', xlabel = P_SDB_GREEN, ylabel = DEPTH)
cal_plot.add_legend(ax)

fig, ax = plt.subplots(1, 1, figsize = (6, 6)) cal_plot = CalibrationPlot(legend_font_size = 15) cal_plot.add_calibration_scatter(calibration, to_cal_p_greens[no_nans], to_cal_depths[no_nans], ax = ax) cal_plot.add_labels(ax, title = f'{P_SDB_GREEN} vs {DEPTH}', xlabel = P_SDB_GREEN, ylabel = DEPTH) cal_plot.add_legend(ax)

Out[13]:

<Axes: title={'center': 'pSDB Green vs DEPTH'}, xlabel='pSDB Green', ylabel='DEPTH'>

SDB Generation¶

In [14]:

SDB_GREEN = 'SDB Green'

multi_image_clip.add_band(SDB_GREEN, calibration.predict(multi_image_clip.select(P_SDB_GREEN)))


fig, axs = plot.get_geofigure(multi_image_clip.crs, 1, 3, figsize = (18, 6))
for image, band, crs, title, vmin, vmax, ax in zip([multi_image_clip, multi_image_clip, in_situ], [P_SDB_GREEN, SDB_GREEN, DEPTH],
                                                   [multi_image_clip.crs, multi_image_clip.crs, in_situ.crs], [P_SDB_GREEN, SDB_GREEN, DEPTH],
                                                   [None, 0, 0], [None, 12, 12], axs):
    ax, mappable = plot.plot_band(image, band, ax, cmap = cmocean.cm.deep, vmin = vmin, vmax = vmax)
    ctx.add_basemap(ax, crs = crs, attribution=False)
    fig.colorbar(mappable, shrink = 0.7)
    plot.add_gridlines(ax)
    ax.set_title(title)

SDB_GREEN = 'SDB Green' multi_image_clip.add_band(SDB_GREEN, calibration.predict(multi_image_clip.select(P_SDB_GREEN))) fig, axs = plot.get_geofigure(multi_image_clip.crs, 1, 3, figsize = (18, 6)) for image, band, crs, title, vmin, vmax, ax in zip([multi_image_clip, multi_image_clip, in_situ], [P_SDB_GREEN, SDB_GREEN, DEPTH], [multi_image_clip.crs, multi_image_clip.crs, in_situ.crs], [P_SDB_GREEN, SDB_GREEN, DEPTH], [None, 0, 0], [None, 12, 12], axs): ax, mappable = plot.plot_band(image, band, ax, cmap = cmocean.cm.deep, vmin = vmin, vmax = vmax) ctx.add_basemap(ax, crs = crs, attribution=False) fig.colorbar(mappable, shrink = 0.7) plot.add_gridlines(ax) ax.set_title(title)

Validation¶

Indexes of values to compare¶

In [15]:

validation_indexes = ~np.isin(np.arange(in_situ.select(DEPTH).size), selection)
validation_indexes

validation_indexes = ~np.isin(np.arange(in_situ.select(DEPTH).size), selection) validation_indexes

Out[15]:

array([ True,  True,  True, ...,  True,  True,  True], shape=(314688,))

Metrics¶

In [16]:

depths = in_situ.select(DEPTH).ravel()[validation_indexes]
green = multi_image_clip.select(SDB_GREEN).ravel()[validation_indexes]

no_nans = is_valid(depths) & is_valid(green)

validation = ValidationSummary(in_situ = depths[no_nans], model = green[no_nans])
validation

depths = in_situ.select(DEPTH).ravel()[validation_indexes] green = multi_image_clip.select(SDB_GREEN).ravel()[validation_indexes] no_nans = is_valid(depths) & is_valid(green) validation = ValidationSummary(in_situ = depths[no_nans], model = green[no_nans]) validation

Out[16]:

N: 133191 | MSD: 0.5760 | MedAE: 0.7906 | Abs_std: 0.99954

Plot¶

In [19]:

fig, axs = plt.subplots(1, 2, figsize = (12, 6))

val_plot = ValidationPlot(legend_font_size = 10)
val_plot.add_densed_scatter(validation, axs[0], x_min = -4, x_max = 30, step = 4, s = 0.1,
                            density = {'method' : 'approximate', 'bins' : 100})
val_plot.add_residuals(validation, axs[1], metrics = ['MedAE', 'Abs_std', 'MSD']) # since we have 100K pixels its better to use approximate gaussian density
axs[0].grid()
axs[1].grid()
fig.tight_layout()

fig, axs = plt.subplots(1, 2, figsize = (12, 6)) val_plot = ValidationPlot(legend_font_size = 10) val_plot.add_densed_scatter(validation, axs[0], x_min = -4, x_max = 30, step = 4, s = 0.1, density = {'method' : 'approximate', 'bins' : 100}) val_plot.add_residuals(validation, axs[1], metrics = ['MedAE', 'Abs_std', 'MSD']) # since we have 100K pixels its better to use approximate gaussian density axs[0].grid() axs[1].grid() fig.tight_layout()

Bathymetry¶

In [18]:

multi_image[SDB_GREEN] = calibration.predict(multi_image[P_SDB_GREEN])
multi_image.mask(multi_image.select(SDB_GREEN) < 10) ## multi_image doens't see more than 10m in this study area so we keep values up-to 10m.

fig, axs = plot.get_geofigure(multi_image.crs, 1, 2, figsize = (18, 6))
axs[0], mappable = plot.plot_band(multi_image, SDB_GREEN, axs[0], cmap = cmocean.cm.deep, vmin = 0, vmax = 12)
ctx.add_basemap(axs[0], crs = multi_image.crs, attribution=False)
fig.colorbar(mappable, label = 'depth (m)')
plot.add_gridlines(axs[0])
axs[0].set_title(SDB_GREEN)

axs[1], mappable = plot.plot_band(in_situ, DEPTH, axs[1], cmap = cmocean.cm.deep, vmin = 0, vmax = 12)
ctx.add_basemap(axs[1], crs = multi_image.crs, attribution=False)
fig.colorbar(mappable, label = 'depth (m)')
plot.add_gridlines(axs[1])
axs[1].set_title('In-Situ')

multi_image[SDB_GREEN] = calibration.predict(multi_image[P_SDB_GREEN]) multi_image.mask(multi_image.select(SDB_GREEN) < 10) ## multi_image doens't see more than 10m in this study area so we keep values up-to 10m. fig, axs = plot.get_geofigure(multi_image.crs, 1, 2, figsize = (18, 6)) axs[0], mappable = plot.plot_band(multi_image, SDB_GREEN, axs[0], cmap = cmocean.cm.deep, vmin = 0, vmax = 12) ctx.add_basemap(axs[0], crs = multi_image.crs, attribution=False) fig.colorbar(mappable, label = 'depth (m)') plot.add_gridlines(axs[0]) axs[0].set_title(SDB_GREEN) axs[1], mappable = plot.plot_band(in_situ, DEPTH, axs[1], cmap = cmocean.cm.deep, vmin = 0, vmax = 12) ctx.add_basemap(axs[1], crs = multi_image.crs, attribution=False) fig.colorbar(mappable, label = 'depth (m)') plot.add_gridlines(axs[1]) axs[1].set_title('In-Situ')

Out[18]:

Text(0.5, 1.0, 'In-Situ')