Satellite Derived Bathymetry - Stumpf Method¶

This example uses stumpf method to derive bathymetry as mentioned in this article https://doi.org/10.4319/lo.2003.48.1_part_2.0547

Imports¶

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
from sensingpy.masks import is_valid
from sensingpy.preprocessing.outliers import IQR
from sensingpy.enums import SENTINEL2_BANDS

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 from sensingpy.masks import is_valid from sensingpy.preprocessing.outliers import IQR from sensingpy.enums import SENTINEL2_BANDS

In [2]:

import warnings

warnings.filterwarnings("ignore")

import warnings warnings.filterwarnings("ignore")

Load Satellite Image¶

In [ ]:

satellite = reader.open(r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181007_S2A_AGV_EXT_29SPA_L2W.nc")
satellite.rename_by_enum(SENTINEL2_BANDS)
satellite.rename({'Rrs_B2' : 'Blue', 'Rrs_B3' : 'Green', 'Rrs_B4' : 'Red'})
satellite.keep_bands( ['Blue', 'Green', 'Red'] )

satellite = reader.open(r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181007_S2A_AGV_EXT_29SPA_L2W.nc") satellite.rename_by_enum(SENTINEL2_BANDS) satellite.rename({'Rrs_B2' : 'Blue', 'Rrs_B3' : 'Green', 'Rrs_B4' : 'Red'}) satellite.keep_bands( ['Blue', 'Green', 'Red'] )

Pseudomodel Green¶

In [4]:

P_SDB_GREEN = 'pSDB Green'

satellite.add_band(P_SDB_GREEN, IQR( stumpf_pseudomodel(satellite.select('Blue'), satellite.select('Green')), distance = 1.5 ))

fig, ax = plot.get_geofigure(satellite.crs, 1, 1)
ax, mappable = plot.plot_band(satellite, P_SDB_GREEN, ax = ax, cmap = cmocean.cm.deep)
ctx.add_basemap(ax, crs = satellite.crs, attribution=False)
fig.colorbar(mappable, label = f'{P_SDB_GREEN} (DII)')
plot.add_gridlines(ax)

P_SDB_GREEN = 'pSDB Green' satellite.add_band(P_SDB_GREEN, IQR( stumpf_pseudomodel(satellite.select('Blue'), satellite.select('Green')), distance = 1.5 )) fig, ax = plot.get_geofigure(satellite.crs, 1, 1) ax, mappable = plot.plot_band(satellite, P_SDB_GREEN, ax = ax, cmap = cmocean.cm.deep) ctx.add_basemap(ax, crs = satellite.crs, attribution=False) fig.colorbar(mappable, label = f'{P_SDB_GREEN} (DII)') plot.add_gridlines(ax)

Out[4]:

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

Load In-Situ Depths¶

In [5]:

DEPTH = 'DEPTH'

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)

DEPTH = 'DEPTH' 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[5]:

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

Align In-Situ with Satellite dimensions¶

In [ ]:

in_situ.align(satellite)
in_situ.dropna()
satellite.clip([in_situ.bbox])

in_situ.align(satellite) in_situ.dropna() satellite.clip([in_situ.bbox])

In-situ and Pseudomodel¶

In [7]:

fig, axs = plot.get_geofigure(satellite.crs, 1, 2)
plot.plot_band(satellite, 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 = satellite.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(satellite.crs, 1, 2) plot.plot_band(satellite, 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 = satellite.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[7]:

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

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 [8]:

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 [9]:

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

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

Fit model¶

In [10]:

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.8242 | y = 55.485x-50.997

Plot model¶

In [11]:

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[11]:

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

SDB Generation¶

In [12]:

SDB_GREEN = 'SDB Green'

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


fig, axs = plot.get_geofigure(satellite.crs, 1, 3, figsize = (18, 6))
for image, band, crs, title, vmin, vmax, ax in zip([satellite, satellite, in_situ], [P_SDB_GREEN, SDB_GREEN, DEPTH],
                                                   [satellite.crs, satellite.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' satellite.add_band(SDB_GREEN, calibration.predict(satellite.select(P_SDB_GREEN))) fig, axs = plot.get_geofigure(satellite.crs, 1, 3, figsize = (18, 6)) for image, band, crs, title, vmin, vmax, ax in zip([satellite, satellite, in_situ], [P_SDB_GREEN, SDB_GREEN, DEPTH], [satellite.crs, satellite.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 [13]:

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[13]:

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

Metrics¶

In [14]:

depths = in_situ.select(DEPTH).ravel()[validation_indexes]
green = satellite.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 = satellite.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[14]:

N: 129496 | MSD: 0.1197 | MedAE: 1.0304 | Abs_std: 1.13966

Plot¶

In [15]:

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()
plt.subplots_adjust(wspace = 0.3)

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() plt.subplots_adjust(wspace = 0.3)

Bathymetry¶

In [16]:

satellite = reader.open(r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181007_S2A_AGV_EXT_29SPA_L2W.nc")
satellite.rename_by_enum(SENTINEL2_BANDS)
satellite.rename({'Rrs_B2' : 'Blue', 'Rrs_B3' : 'Green', 'Rrs_B4' : 'Red'})
satellite.keep_bands(['Blue', 'Green', 'Red'])

satellite.add_band(P_SDB_GREEN, IQR( stumpf_pseudomodel(satellite.select('Blue'), satellite.select('Green')), distance = 1.5 ))
satellite.add_band(SDB_GREEN, calibration.predict(satellite.select(P_SDB_GREEN)))
satellite.mask(satellite.select(SDB_GREEN) < 10) ## Satellite doens't see more than 10m in this study area so we keep values up-to 10m.


fig, axs = plot.get_geofigure(satellite.crs, 1, 2, figsize = (18, 6))

axs[0], mappable = plot.plot_band(satellite, SDB_GREEN, axs[0], cmap = cmocean.cm.deep, vmin = 0, vmax = 12)
ctx.add_basemap(axs[0], crs = satellite.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 = satellite.crs, attribution=False)
fig.colorbar(mappable, label = 'depth (m)')
plot.add_gridlines(axs[1])
axs[1].set_title('In-Situ')

satellite = reader.open(r"E:\batimetria\algarve\satelite\local_original\Pleamar\20181007_S2A_AGV_EXT_29SPA_L2W.nc") satellite.rename_by_enum(SENTINEL2_BANDS) satellite.rename({'Rrs_B2' : 'Blue', 'Rrs_B3' : 'Green', 'Rrs_B4' : 'Red'}) satellite.keep_bands(['Blue', 'Green', 'Red']) satellite.add_band(P_SDB_GREEN, IQR( stumpf_pseudomodel(satellite.select('Blue'), satellite.select('Green')), distance = 1.5 )) satellite.add_band(SDB_GREEN, calibration.predict(satellite.select(P_SDB_GREEN))) satellite.mask(satellite.select(SDB_GREEN) < 10) ## Satellite doens't see more than 10m in this study area so we keep values up-to 10m. fig, axs = plot.get_geofigure(satellite.crs, 1, 2, figsize = (18, 6)) axs[0], mappable = plot.plot_band(satellite, SDB_GREEN, axs[0], cmap = cmocean.cm.deep, vmin = 0, vmax = 12) ctx.add_basemap(axs[0], crs = satellite.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 = satellite.crs, attribution=False) fig.colorbar(mappable, label = 'depth (m)') plot.add_gridlines(axs[1]) axs[1].set_title('In-Situ')

Out[16]:

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