Analysing snow cover in the Alps using Sentinel-2

from datetime import datetime
import xarray as xr
import dask
from pyproj import Transformer
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
import numpy as np
import pystac_client

catchment_outline = gpd.read_file("./geojson/catchment_outline_cc_31.geojson")
bbox = list(catchment_outline.total_bounds)  # minx, miny, maxx, maxy
bbox

catalog = pystac_client.Client.open("https://stac.core.eopf.eodc.eu")
search = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=bbox,
    datetime="2018-02-01/2018-06-30",
    query={"eo:cloud_cover": {"lt": 75}},
)

items = list(search.items())
print(f"{len(items)} STAC Items found!")

def reproject_bbox(bbox, src_crs="EPSG:4326", dst_crs="EPSG:32632"):
    transformer = Transformer.from_crs(src_crs, dst_crs, always_xy=True)
    xmin, ymin = transformer.transform(bbox[0], bbox[1])
    xmax, ymax = transformer.transform(bbox[2], bbox[3])
    return [xmin, ymin, xmax, ymax]


def load_tile(item, bbox_ll):
    href = item.assets["product"].href
    ds = xr.open_datatree(href, engine="zarr", chunks={}, mask_and_scale=True)

    dst_crs = item.properties.get("proj:code")
    if dst_crs is None:
        raise ValueError(f"No proj:code in item {item.id}")

    utm_bbox = reproject_bbox(bbox_ll, dst_crs=dst_crs)

    ds_ref = ds["measurements"]["reflectance"]["r20m"].to_dataset()[
        ["b02", "b03", "b04", "b11"]
    ]
    ds_scl = ds["conditions"]["mask"]["l2a_classification"]["r20m"].to_dataset()
    ds_out = (
        xr.merge([ds_ref, ds_scl])
        .sel(x=slice(utm_bbox[0], utm_bbox[2]), y=slice(utm_bbox[3], utm_bbox[1]))
        .expand_dims(time=[item.datetime])
    )
    return ds_out


delayed_results = [dask.delayed(load_tile)(item, bbox) for item in items]
results = dask.compute(*delayed_results)

# Get CRS from STAC item
crs = items[0].properties.get("proj:code")
datacube = xr.concat(results, dim="time")
datacube = datacube.assign_coords(time=pd.to_datetime(datacube.time.values)).sortby(
    "time"
)
datacube = datacube.groupby(datacube.time.dt.floor("D")).max()
datacube = datacube.rename(floor="time").sortby("time")

datacube = datacube.rio.write_crs(crs, inplace=True)
datacube

datacube = datacube.sortby("y", ascending=False)
datacube = datacube.rio.write_transform()
catchment_outline = catchment_outline.to_crs(datacube.rio.crs)

datacube = datacube.rio.clip(
    catchment_outline.geometry.values,
    catchment_outline.crs,
    all_touched=True,
    drop=True,
)
no_data_mask = ~np.isnan(datacube.b02)
datacube

datacube_sample = datacube.isel(time=15)
rgb = datacube_sample[["b04", "b03", "b02"]].to_dataarray(dim="band")
rgb = rgb.assign_coords(band=["red", "green", "blue"])

vmin = float(rgb.min())
vmax = float(rgb.quantile(0.95))
rgb_scaled = (rgb - vmin) / (vmax - vmin)
rgb_scaled = rgb_scaled.clip(0, 1)

rgb_img = rgb_scaled.transpose("y", "x", "band").values

rgb_img = np.ma.masked_invalid(rgb_img)

x = datacube_sample.x.values
y = datacube_sample.y.values

fig, ax = plt.subplots(figsize=(8, 8))
ax.imshow(rgb_img, extent=[x.min(), x.max(), y.min(), y.max()])

catchment_outline.boundary.plot(ax=ax, edgecolor="yellow", linewidth=2)
ax.set_title("RGB composite with catchment boundary overlay")
ax.axis("off")
plt.show()

ndsi = (datacube["b03"] - datacube["b11"]) / (datacube["b03"] + datacube["b11"])
ndsi = ndsi.rename("ndsi")
ndsi

ndsi_mask = ndsi > 0.42
snowmap = ndsi_mask.rename("snow_map").where(no_data_mask)
snowmap

cloud_mask = datacube["scl"].isin([3, 8, 9])
snowmap_cloudfree = snowmap.where(~cloud_mask, other=2)
snowmap_cloudfree

snowmap_1d = snowmap_cloudfree.sel(time=slice("2018-02-10", "2018-02-12"))
fig, ax = plt.subplots(figsize=(10, 8))
snowmap_1d.isel(time=0).plot(ax=ax, cmap="viridis")
ax.set_title(f"Snow Map - {snowmap_1d.time.values[0]}")
plt.show()

snow_pixels = (snowmap_cloudfree == 1).where(no_data_mask)
cloud_pixels = cloud_mask.where(no_data_mask)
valid_pixels = (snowmap_cloudfree >= 0).where(no_data_mask)

pixel_maps = xr.Dataset(
    {
        "snow_pixels": snow_pixels,
        "cloud_pixels": cloud_pixels,
        "valid_pixels": valid_pixels,
    }
)

catchment_stats = pixel_maps.sum(dim=["x", "y"])

print("Catchment statistics over time:")
print(catchment_stats)

total_pixels = catchment_stats["valid_pixels"]
snow_percentage = catchment_stats["snow_pixels"] / total_pixels * 100
cloud_percentage = catchment_stats["cloud_pixels"] / total_pixels * 100

catchment_stats = catchment_stats.assign(
    snow_percentage=snow_percentage, cloud_percentage=cloud_percentage
).compute()

cloud_threshold = 25

cloud_percent = catchment_stats["cloud_percentage"]

valid_mask = cloud_percent <= cloud_threshold

valid_timesteps = catchment_stats.where(valid_mask, drop=True)

snowmap_filtered = snowmap_cloudfree.where(valid_mask, drop=True)

print("Original timesteps:", len(catchment_stats.time))
print(f"Timesteps with <{cloud_threshold}% clouds:", len(valid_timesteps.time))

plt.figure(figsize=(12, 5))

cloud_percent = cloud_percent.assign_coords(
    time=pd.to_datetime(cloud_percent.time.astype("int64"))
)

cloud_percent.plot(marker="o", label="Cloud Percentage")

plt.axhline(y=25, color="r", linestyle="--", label="25% threshold")

plt.title("Cloud Percentage Over Time with Threshold")
plt.ylabel("Cloud Percentage (%)")
plt.xlabel("Time")
plt.legend()
plt.grid(True)

plt.show()

snow_percentage = catchment_stats["snow_percentage"]

plt.figure(figsize=(12, 5))

snow_percentage = snow_percentage.assign_coords(
    time=pd.to_datetime(snow_percentage.time.astype("int64"))
)

snow_percentage.plot(marker="o", label="Snow Percentage")

plt.title("Snow Percentage Over Time with Threshold")
plt.ylabel("Snow Percentage (%)")
plt.xlabel("Time")
plt.legend()
plt.grid(True)

plt.show()

filtered_snow_percentage = snow_percentage.where(valid_mask, drop=True)

filtered_snow_percentage = filtered_snow_percentage.assign_coords(
    time=pd.to_datetime(filtered_snow_percentage.time.astype("int64"))
)

plt.figure(figsize=(12, 5))

filtered_snow_percentage.plot(marker="o", label="Snow Percentage")

plt.title("Cloud-Filtered Snow Percentage Over Time")
plt.ylabel("Snow Percentage (%)")
plt.xlabel("Time")
plt.legend()
plt.grid(True)
plt.show()

discharge_ds = pd.read_csv(
    "./csv/ADO_DSC_ITH1_0025_CC_31.csv", sep=",", index_col="Time", parse_dates=True
)
discharge_ds.head()

start_date = datetime(2018, 2, 1)
end_date = datetime(2018, 6, 30)

snow_percent = filtered_snow_percentage.sel(time=slice(start_date, end_date))
discharge_ds = discharge_ds.loc[start_date:end_date]

fig, ax0 = plt.subplots(1, figsize=(10, 5), sharey=True)
ax1 = ax0.twinx()

discharge_ds.discharge_m3_s.plot(
    label="Discharge", ylabel="Discharge (m$^3$/s)", ax=ax0
)
snow_percent.plot(marker="o", ax=ax1, color="orange")

ax0.legend(loc="center left", bbox_to_anchor=(0, 0.6))
ax1.set_ylabel("Snow cover area (%)")
ax1.legend(loc="center left", bbox_to_anchor=(0, 0.5), labels=["SCA"])

plt.show()