Remote Sensing for Water Bodies¶

Imports¶

In [1]:

import matplotlib.pyplot as plt
import geopandas as gpd
import numpy as np


from sensingpy import reader, plot
from sensingpy.enums import SENTINEL2_BANDS
from sensingpy.preprocessing.outliers import IQR

import matplotlib.pyplot as plt import geopandas as gpd import numpy as np from sensingpy import reader, plot from sensingpy.enums import SENTINEL2_BANDS from sensingpy.preprocessing.outliers import IQR

In [2]:

import warnings

warnings.filterwarnings("ignore")

import warnings warnings.filterwarnings("ignore")

Load Satellite Image¶

In [ ]:

satellite = reader.open(r"D:\marine_observatory\raw\Sentinel-2\30SXG\S2C_MSI_2025_03_25_11_01_08_T30SXG_L2R.nc")
satellite.rename_by_enum(SENTINEL2_BANDS)

satellite = reader.open(r"D:\marine_observatory\raw\Sentinel-2\30SXG\S2C_MSI_2025_03_25_11_01_08_T30SXG_L2R.nc") satellite.rename_by_enum(SENTINEL2_BANDS)

In [ ]:

mmm = gpd.read_file(r'D:\batimetria\scripts\areas\mmn.json').to_crs(satellite.crs).geometry

satellite.clip(mmm)

mmm = gpd.read_file(r'D:\batimetria\scripts\areas\mmn.json').to_crs(satellite.crs).geometry satellite.clip(mmm)

RGB¶

In [5]:

fig, ax = plot.get_geofigure(satellite.crs, 1, 1)
ax = plot.plot_rgb(satellite, 'rhos_B4', 'rhos_B3', 'rhos_B2', ax = ax, brightness = 5)

fig, ax = plot.get_geofigure(satellite.crs, 1, 1) ax = plot.plot_rgb(satellite, 'rhos_B4', 'rhos_B3', 'rhos_B2', ax = ax, brightness = 5)

False Color (Vegetation)¶

In [6]:

fig, ax = plot.get_geofigure(satellite.crs, 1, 1)
ax = plot.plot_rgb(satellite, 'rhos_B8', 'rhos_B4', 'rhos_B3', ax = ax, brightness = 5)

fig, ax = plot.get_geofigure(satellite.crs, 1, 1) ax = plot.plot_rgb(satellite, 'rhos_B8', 'rhos_B4', 'rhos_B3', ax = ax, brightness = 5)

Water Index (NDWI)¶

In [ ]:

satellite.add_band('NDWI', IQR(satellite.normalized_diference('rhos_B3', 'rhos_B8')))


fig, ax = plot.get_geofigure(satellite.crs, 1, 1)
ax, mappable = plot.plot_band(satellite, 'NDWI', ax = ax, vmin = -1, vmax = 1, cmap = 'GnBu')
fig.colorbar(mappable)

satellite.add_band('NDWI', IQR(satellite.normalized_diference('rhos_B3', 'rhos_B8'))) fig, ax = plot.get_geofigure(satellite.crs, 1, 1) ax, mappable = plot.plot_band(satellite, 'NDWI', ax = ax, vmin = -1, vmax = 1, cmap = 'GnBu') fig.colorbar(mappable)

Out[ ]:

<matplotlib.colorbar.Colorbar at 0x22314f5a5d0>

Mask Land using SWIR 1 band¶

In [ ]:

satellite.mask(satellite.select('rhot_B11') < 0.1)

satellite.mask(satellite.select('rhot_B11') < 0.1)

Transform Water Surface Reflectance to Rrs (Water Leaving)¶

In [10]:

for band in [band for band in satellite.band_names if 'rhos' in band]:
    satellite.add_band(band.replace('rhos', 'Rrs'), satellite.select(band) / np.pi)

for band in [band for band in satellite.band_names if 'rhos' in band]: satellite.add_band(band.replace('rhos', 'Rrs'), satellite.select(band) / np.pi)

Spectral analysis¶

In [17]:

points = gpd.read_file(r"D:\batimetria\scripts\areas\mmn_points\mmn_points.shp").to_crs(satellite.crs).geometry.iloc[4::2]


fig, ax = plot.get_geofigure(satellite.crs, 1, 1, figsize = (24, 6))
ax = plot.plot_rgb(satellite, 'rhos_B4', 'rhos_B3', 'rhos_B2', ax = ax, brightness = 5)

for idx, point in enumerate(points):
    ax.scatter(point.x, point.y, s = 50, label = f'P{idx + 1}')
    ax.text(point.x, point.y + 1000, f'P{idx + 1}', fontsize = 12, color = 'white')

points = gpd.read_file(r"D:\batimetria\scripts\areas\mmn_points\mmn_points.shp").to_crs(satellite.crs).geometry.iloc[4::2] fig, ax = plot.get_geofigure(satellite.crs, 1, 1, figsize = (24, 6)) ax = plot.plot_rgb(satellite, 'rhos_B4', 'rhos_B3', 'rhos_B2', ax = ax, brightness = 5) for idx, point in enumerate(points): ax.scatter(point.x, point.y, s = 50, label = f'P{idx + 1}') ax.text(point.x, point.y + 1000, f'P{idx + 1}', fontsize = 12, color = 'white')

In [ ]:

bands = [band for band in satellite.band_names if 'Rrs' in band]
spectral = satellite.extract_values(points.x.values, points.y.values, bands = bands)

fig, ax = plt.subplots(1, 1, figsize = (12, 6))
for i in range(len(points)):
    ax.plot(bands, spectral[:, i], '-o', label = f'P{i + 1}')
ax.set_xlabel('Wavelength (nm)')
ax.set_ylabel('Reflectance')
ax.legend()
ax.grid()

bands = [band for band in satellite.band_names if 'Rrs' in band] spectral = satellite.extract_values(points.x.values, points.y.values, bands = bands) fig, ax = plt.subplots(1, 1, figsize = (12, 6)) for i in range(len(points)): ax.plot(bands, spectral[:, i], '-o', label = f'P{i + 1}') ax.set_xlabel('Wavelength (nm)') ax.set_ylabel('Reflectance') ax.legend() ax.grid()