Dielectric Properties

import json

import holoviews as hv
import intake
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import numpy as np
import rioxarray  # noqa: F401
import xarray as xr
from holoviews.streams import RangeXY
from matplotlib.colors import BoundaryNorm, ListedColormap

hv.extension("bokeh")

uri = "https://git.geo.tuwien.ac.at/public_projects/microwave-remote-sensing/-/raw/main/microwave-remote-sensing.yml"
cat = intake.open_catalog(uri)
sig0_da = cat.neusiedler.read().sig0.compute()

sig0_da.isel(time=0).plot(robust=True, cmap="Greys_r").axes.set_aspect("equal")

cor_da = cat.corine.read().land_cover.compute()

# Load encoding
with cat.corine_cmap.read()[0] as f:
    color_mapping_data = json.load(f)

# Get mapping
color_mapping = {item["value"]: item for item in color_mapping_data["land_cover"]}

# Create cmap and norm for plotting
colors = [info["color"] for info in color_mapping.values()]
categories = [info["value"] for info in color_mapping.values()]
cmap = ListedColormap(colors)
norm = BoundaryNorm(categories + [max(categories) + 1], len(categories))

# Get landcover codes present in the image
present_landcover_codes = np.unique(cor_da.values[~np.isnan(cor_da.values)].astype(int))

# Get colors + text for legend
handles = [
    mpatches.Patch(color=info["color"], label=(f"{info['value']} - " + (info["label"])))
    for info in color_mapping.values()
    if info["value"] in present_landcover_codes
]

# Create the plot
cor_da.plot(figsize=(10, 10), cmap=cmap, norm=norm, add_colorbar=False).axes.set_aspect(
    "equal"
)

plt.legend(
    handles=handles,
    bbox_to_anchor=(1.01, 1),
    loc="upper left",
    borderaxespad=0,
    fontsize=7,
)
plt.title("CORINE Land Cover (EPSG:27704)")

var_ds = xr.merge([sig0_da, cor_da]).drop_vars("band")
var_ds

# 41 = encoded value for water bodies
waterbodies_mask = var_ds.land_cover == 41
waterbodies_mask.plot().axes.set_aspect("equal")

waterbodies_sig0 = var_ds.sig0.isel(time=0).where(waterbodies_mask)
waterbodies_sig0.plot(robust=True, cmap="Greys_r").axes.set_aspect("equal")

waterbodies_sig0.plot.hist(bins=50, edgecolor="black")

robust_min = var_ds.sig0.quantile(0.02).item()
robust_max = var_ds.sig0.quantile(0.98).item()

bin_edges = [
    i + j * 0.5
    for i in range(int(robust_min) - 2, int(robust_max) + 2)
    for j in range(2)
]

land_cover = {"\xa0\xa0\xa0 Complete Land Cover": 1}
land_cover.update(
    {
        f"{int(value): 02} {color_mapping[value]['label']}": int(value)
        for value in present_landcover_codes
    }
)
time = var_ds.sig0["time"].values

rangexy = RangeXY()


def load_image(time, land_cover, x_range, y_range):
    """Callback Function Landcover.

    Parameters
    ----------
    time: panda.datatime
        time slice
    landcover: int
        land cover type
    x_range: array_like
        longitude range
    y_range: array_like
        latitude range

    Returns
    -------
    holoviews.Image

    """
    if land_cover == "\xa0\xa0\xa0 Complete Land Cover":
        sig0_selected_ds = var_ds.sig0.sel(time=time)

    else:
        land_cover_value = int(land_cover.split()[0])
        mask_ds = var_ds.land_cover == land_cover_value
        sig0_selected_ds = var_ds.sig0.sel(time=time).where(mask_ds)

    hv_ds = hv.Dataset(sig0_selected_ds)
    img = hv_ds.to(hv.Image, ["x", "y"])

    if x_range and y_range:
        img = img.select(x=x_range, y=y_range)

    return hv.Image(img)


dmap = (
    hv.DynamicMap(load_image, kdims=["Time", "Landcover"], streams=[rangexy])
    .redim.values(Time=time, Landcover=land_cover)
    .hist(normed=True, bins=bin_edges)
)

image_opts = hv.opts.Image(
    cmap="Greys_r",
    colorbar=True,
    tools=["hover"],
    clim=(robust_min, robust_max),
    aspect="equal",
    framewise=False,
    frame_height=500,
    frame_width=500,
)

hist_opts = hv.opts.Histogram(width=350, height=555)

dmap.opts(image_opts, hist_opts)