{
"cells": [
{
"cell_type": "markdown",
"id": "0",
"metadata": {},
"source": [
"# Dielectric Properties\n",
"\n",
"\n",
"\n",
"In this notebook, we will investigate the varying backscatter values associated with different land surfaces like water bodies, forests, grasslands and urban areas. We will use backscatter data from the Sentinel-1 satellite and we will utilize the CORINE Land Cover dataset to classify and extrapolate these surfaces, enabling us to analyze how different land cover types influence backscatter responses.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1",
"metadata": {},
"outputs": [],
"source": [
"import json\n",
"\n",
"import holoviews as hv\n",
"import intake\n",
"import matplotlib.patches as mpatches\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import rioxarray # noqa: F401\n",
"import xarray as xr\n",
"from holoviews.streams import RangeXY\n",
"from matplotlib.colors import BoundaryNorm, ListedColormap\n",
"\n",
"hv.extension(\"bokeh\")"
]
},
{
"cell_type": "markdown",
"id": "2",
"metadata": {},
"source": [
"## Load Sentinel-1 Data\n",
"\n",
"For our analysis we are using sigma naught backscatering data from Sentinel-1. The images we are analyzing cover the region south of Vienna and west of Lake Neusiedl. We load the data and and apply again a preprocessing function. Here we extract the scaling factor and the date the image was taken from the metadata. We will focus our attention to a smaller area containing a part of the Lake Neusiedl Lake and its surrounding land. The obtained`xarray` dataset and is then converted to an array, because we only have one variable, the VV backscatter values.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3",
"metadata": {},
"outputs": [],
"source": [
"url = \"https://huggingface.co/datasets/martinschobben/microwave-remote-sensing/resolve/main/microwave-remote-sensing.yml\"\n",
"cat = intake.open_catalog(url)\n",
"sig0_da = cat.neusiedler.read().sig0.compute()"
]
},
{
"cell_type": "markdown",
"id": "4",
"metadata": {},
"source": [
"Let's have a look at the data by plotting the first timeslice.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5",
"metadata": {},
"outputs": [],
"source": [
"sig0_da.isel(time=0).plot(robust=True, cmap=\"Greys_r\").axes.set_aspect(\"equal\")"
]
},
{
"cell_type": "markdown",
"id": "6",
"metadata": {},
"source": [
"## Load CORINE Landcover Data\n",
"\n",
"We will load the CORINE Land Cover, which is a pan-European land cover and land use inventory with 44 thematic classes. The resolution of this classification is 100 by 100m and the file was created in 2018\n",
"([CORINE Land Cover](https://land.copernicus.eu/en/products/corine-land-cover)).\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7",
"metadata": {},
"outputs": [],
"source": [
"cor_da = cat.corine.read().land_cover.compute()"
]
},
{
"cell_type": "markdown",
"id": "8",
"metadata": {},
"source": [
"### Colormapping and Encoding\n",
"\n",
"For the different land cover types we use the official color encoding.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9",
"metadata": {},
"outputs": [],
"source": [
"# Load encoding\n",
"with cat.corine_cmap.read()[0] as f:\n",
" color_mapping_data = json.load(f)\n",
"\n",
"# Get mapping\n",
"color_mapping = {item[\"value\"]: item for item in color_mapping_data[\"land_cover\"]}\n",
"\n",
"# Create cmap and norm for plotting\n",
"colors = [info[\"color\"] for info in color_mapping.values()]\n",
"categories = [info[\"value\"] for info in color_mapping.values()]\n",
"cmap = ListedColormap(colors)\n",
"norm = BoundaryNorm(categories + [max(categories) + 1], len(categories))"
]
},
{
"cell_type": "markdown",
"id": "10",
"metadata": {},
"source": [
"Now we can plot the CORINE Land Cover dataset.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "11",
"metadata": {},
"outputs": [],
"source": [
"# Get landcover codes present in the image\n",
"present_landcover_codes = np.unique(cor_da.values[~np.isnan(cor_da.values)].astype(int))\n",
"\n",
"# Get colors + text for legend\n",
"handles = [\n",
" mpatches.Patch(color=info[\"color\"], label=(f'{info[\"value\"]} - ' + (info[\"label\"])))\n",
" for info in color_mapping.values()\n",
" if info[\"value\"] in present_landcover_codes\n",
"]\n",
"\n",
"# Create the plot\n",
"cor_da.plot(figsize=(10, 10), cmap=cmap, norm=norm, add_colorbar=False).axes.set_aspect(\n",
" \"equal\"\n",
")\n",
"\n",
"plt.legend(\n",
" handles=handles,\n",
" bbox_to_anchor=(1.01, 1),\n",
" loc=\"upper left\",\n",
" borderaxespad=0,\n",
" fontsize=7,\n",
")\n",
"plt.title(\"CORINE Land Cover (EPSG:27704)\")"
]
},
{
"cell_type": "markdown",
"id": "12",
"metadata": {},
"source": [
"Now we are ready to merge the backscatter data (`sig0_da`) with the land cover dataset (`cor_da`) to have one dataset combining all data.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "13",
"metadata": {},
"outputs": [],
"source": [
"var_ds = xr.merge([sig0_da, cor_da]).drop_vars(\"band\")\n",
"var_ds"
]
},
{
"cell_type": "markdown",
"id": "14",
"metadata": {},
"source": [
"## Backscatter Variability\n",
"\n",
"With this combined dataset we can study backscatter variability in relation to natural media. For example we can look at the backscatter variability for water by clipping the dataset to only contain the land cover class water, like so:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "15",
"metadata": {},
"outputs": [],
"source": [
"# 41 = encoded value for water bodies\n",
"waterbodies_mask = var_ds.land_cover == 41\n",
"waterbodies_mask.plot().axes.set_aspect(\"equal\")"
]
},
{
"cell_type": "markdown",
"id": "16",
"metadata": {},
"source": [
"This gives use backscatter values over water only.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "17",
"metadata": {},
"outputs": [],
"source": [
"waterbodies_sig0 = var_ds.sig0.isel(time=0).where(waterbodies_mask)\n",
"waterbodies_sig0.plot(robust=True, cmap=\"Greys_r\").axes.set_aspect(\"equal\")"
]
},
{
"cell_type": "markdown",
"id": "18",
"metadata": {},
"source": [
"To get an idea of the variability we can create a histogram. Radar backscatter from water bodies fluctuates with surface roughness, which changes with wind conditions, creating spatial and temporal variations in signal intensity.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "19",
"metadata": {},
"outputs": [],
"source": [
"waterbodies_sig0.plot.hist(bins=50, edgecolor=\"black\")"
]
},
{
"cell_type": "markdown",
"id": "20",
"metadata": {},
"source": [
"## Variability over Time\n",
"\n",
"Next we will look at the changes in variability in backscatter values over time for each of the CORINE Land Cover types. We do this by creating the following interactive plot. We can spot that backscatter in agricultural fields varies due to seasonal cycles like planting, growing, and harvesting, each of which changes vegetation structure. Changes in backscatter are strongly related to soil moisture content from irrigation or rainfall. Ultimately, phenological stages of crops and canopy moisture dynamics can affect the backscatter signal.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "21",
"metadata": {},
"outputs": [],
"source": [
"robust_min = var_ds.sig0.quantile(0.02).item()\n",
"robust_max = var_ds.sig0.quantile(0.98).item()\n",
"\n",
"bin_edges = [\n",
" i + j * 0.5\n",
" for i in range(int(robust_min) - 2, int(robust_max) + 2)\n",
" for j in range(2)\n",
"]\n",
"\n",
"land_cover = {\"\\xa0\\xa0\\xa0 Complete Land Cover\": 1}\n",
"land_cover.update(\n",
" {\n",
" f\"{int(value): 02} {color_mapping[value]['label']}\": int(value)\n",
" for value in present_landcover_codes\n",
" }\n",
")\n",
"time = var_ds.sig0[\"time\"].values\n",
"\n",
"rangexy = RangeXY()\n",
"\n",
"\n",
"def load_image(time, land_cover, x_range, y_range):\n",
" \"\"\"\n",
" Callback Function Landcover.\n",
"\n",
" Parameters\n",
" ----------\n",
" time: panda.datatime\n",
" time slice\n",
" landcover: int\n",
" land cover type\n",
" x_range: array_like\n",
" longitude range\n",
" y_range: array_like\n",
" latitude range\n",
"\n",
" Returns\n",
" -------\n",
" holoviews.Image\n",
" \"\"\"\n",
"\n",
" if land_cover == \"\\xa0\\xa0\\xa0 Complete Land Cover\":\n",
" sig0_selected_ds = var_ds.sig0.sel(time=time)\n",
"\n",
" else:\n",
" land_cover_value = int(land_cover.split()[0])\n",
" mask_ds = var_ds.land_cover == land_cover_value\n",
" sig0_selected_ds = var_ds.sig0.sel(time=time).where(mask_ds)\n",
"\n",
" hv_ds = hv.Dataset(sig0_selected_ds)\n",
" img = hv_ds.to(hv.Image, [\"x\", \"y\"])\n",
"\n",
" if x_range and y_range:\n",
" img = img.select(x=x_range, y=y_range)\n",
"\n",
" return hv.Image(img)\n",
"\n",
"\n",
"dmap = (\n",
" hv.DynamicMap(load_image, kdims=[\"Time\", \"Landcover\"], streams=[rangexy])\n",
" .redim.values(Time=time, Landcover=land_cover)\n",
" .hist(normed=True, bins=bin_edges)\n",
")\n",
"\n",
"image_opts = hv.opts.Image(\n",
" cmap=\"Greys_r\",\n",
" colorbar=True,\n",
" tools=[\"hover\"],\n",
" clim=(robust_min, robust_max),\n",
" aspect=\"equal\",\n",
" framewise=False,\n",
" frame_height=500,\n",
" frame_width=500,\n",
")\n",
"\n",
"hist_opts = hv.opts.Histogram(width=350, height=555)\n",
"\n",
"dmap.opts(image_opts, hist_opts)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}