Spectral Clustering

Overview

The current notebook will demonstrate a simplified machine learning approach to observe the change in a lake water’s extent across time. In order to identify the water, we can use spectral clustering to classify each grid cell into a category based on the similarity of the combined set of pixels across wavelength-bands in our image stacks.

Our example approach uses a version of spectral clustering from dask_ml that is a scalable equivalent of what is available in scikit-learn. We will begin this approach with a single image stack and then conduct a direct comparison on the results from different time points.

This workflow uses data from Microsoft Planetary Computer but it can be adapted to work with any data ingestion approach from this cookbook.

Prerequisites

Concepts	Importance	Notes
Data Ingestion - Geospatial-Specific Tooling	Necessary
scikit-learn	Helpful	Spectral clustering
dask_ml	Helpful	Spectral clustering at scale

Time to learn: 20 minutes.

Imports

import numpy as np
import odc.stac
import pandas as pd
import planetary_computer
import pystac_client
import xarray as xr
from dask.distributed import Client
from pystac.extensions.eo import EOExtension as eo
from dask_ml.cluster import SpectralClustering
import pyproj

# Viz
import hvplot.xarray

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/dask/dataframe/_pyarrow_compat.py:17: FutureWarning: Minimal version of pyarrow will soon be increased to 14.0.1. You are using 12.0.1. Please consider upgrading.
  warnings.warn(

Loading Data

Let’s start by loading some Landsat data. These steps are covered in the Data Ingestion - Planetary Computer prerequisite.

Search the catalog

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

bbox = [-118.89, 38.54, -118.57, 38.84]  # Region over a lake in Nevada, USA
datetime = "2017-06-01/2017-09-30"  # Summer months of 2017
collection = "landsat-c2-l2"
platform = "landsat-8"
cloudy_less_than = 1  # percent

search = catalog.search(
    collections=["landsat-c2-l2"],
    bbox=bbox,
    datetime=datetime,
    query={"eo:cloud_cover": {"lt": cloudy_less_than}, "platform": {"in": [platform]}},
)
items = search.get_all_items()
print(f"Returned {len(items)} Items:")
[[i, item.id] for i, item in enumerate(items)]

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/pystac_client/item_search.py:849: FutureWarning: get_all_items() is deprecated, use item_collection() instead.
  warnings.warn(

Returned 3 Items:

[[0, 'LC08_L2SP_042033_20170718_02_T1'],
 [1, 'LC08_L2SP_042033_20170702_02_T1'],
 [2, 'LC08_L2SP_042033_20170616_02_T1']]

Load a dataset

item = items[1]  # select one of the results

assets = []
for _, asset in item.assets.items():
    try:
        assets.append(asset.extra_fields["eo:bands"][0])
    except:
        pass

cols_ordered = [
    "common_name",
    "description",
    "name",
    "center_wavelength",
    "full_width_half_max",
]
bands = pd.DataFrame.from_dict(assets)[cols_ordered]
bands

	common_name	description	name	center_wavelength	full_width_half_max
0	red	Visible red	OLI_B4	0.65	0.04
1	blue	Visible blue	OLI_B2	0.48	0.06
2	green	Visible green	OLI_B3	0.56	0.06
3	nir08	Near infrared	OLI_B5	0.87	0.03
4	lwir11	Long-wave infrared	TIRS_B10	10.90	0.59
5	swir16	Short-wave infrared	OLI_B6	1.61	0.09
6	swir22	Short-wave infrared	OLI_B7	2.20	0.19
7	coastal	Coastal/Aerosol	OLI_B1	0.44	0.02

ds_2017 = odc.stac.stac_load(
    [item],
    bands=bands.common_name.values,
    bbox=bbox,
    chunks={},  # <-- use Dask
).isel(time=0)

Retain CRS Attribute

epsg = item.properties["proj:epsg"]
ds_2017.attrs["crs"] = f"epsg:{epsg}"

da_2017 = ds_2017.to_array(dim="band")
da_2017

<xarray.DataArray (band: 8, y: 1128, x: 950)>
dask.array<stack, shape=(8, 1128, 950), dtype=uint16, chunksize=(1, 1128, 950), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
  * x            (x) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.638e+05
    spatial_ref  int32 32611
    time         datetime64[ns] 2017-07-02T18:33:06.200763
  * band         (band) object 'red' 'blue' 'green' ... 'swir22' 'coastal'
Attributes:
    crs:      epsg:32611

Reshaping Data

The shape of our data is currently n_bands, n_y, n_x. In order for dask-ml / scikit-learn to consume our data, we’ll need to reshape our image stacks into n_samples, n_features, where n_features is the number of wavelength-bands and n_samples is the total number of pixels in each wavelength-band image. Essentially, we’ll be creating a vector of pixels out of each image, where each pixel has multiple features (bands), but the ordering of the pixels is no longer relevant to the computation.

By using xarray methods to flatten the data, we can keep track of the coordinate labels ‘x’ and ‘y’ along the way. This means that we have the ability to reshape back to our original array at any time with no information loss!

flattened_xda = da_2017.stack(z=("x", "y"))  # flatten each band
flattened_t_xda = flattened_xda.transpose("z", "band")
flattened_t_xda

<xarray.DataArray (z: 1071600, band: 8)>
dask.array<transpose, shape=(1071600, 8), dtype=uint16, chunksize=(1071600, 1), chunktype=numpy.ndarray>
Coordinates:
    spatial_ref  int32 32611
    time         datetime64[ns] 2017-07-02T18:33:06.200763
  * band         (band) object 'red' 'blue' 'green' ... 'swir22' 'coastal'
  * z            (z) object MultiIndex
  * x            (z) float64 3.353e+05 3.353e+05 ... 3.638e+05 3.638e+05
  * y            (z) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
Attributes:
    crs:      epsg:32611

Standardize Data

Now that we have the data in the correct shape, let’s standardize (or rescale) the values of the data. We do this to get all the flattened image vectors onto a common scale while preserving the differences in the ranges of values. Again, we’ll demonstrate doing this first in NumPy and then xarray.

with xr.set_options(keep_attrs=True):
    rescaled_xda = (flattened_t_xda - flattened_t_xda.mean()) / flattened_t_xda.std()
rescaled_xda

<xarray.DataArray (z: 1071600, band: 8)>
dask.array<truediv, shape=(1071600, 8), dtype=float64, chunksize=(1071600, 1), chunktype=numpy.ndarray>
Coordinates:
    spatial_ref  int32 32611
    time         datetime64[ns] 2017-07-02T18:33:06.200763
  * band         (band) object 'red' 'blue' 'green' ... 'swir22' 'coastal'
  * z            (z) object MultiIndex
  * x            (z) float64 3.353e+05 3.353e+05 ... 3.638e+05 3.638e+05
  * y            (z) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
Attributes:
    crs:      epsg:32611

Info

Above, we are using a context manager “with xr.set_options(keep_attrs=True):” to retain the array’s attributes through the operations. That is, we want any metadata like ‘crs’ to stay with our result so we can use ‘geo=True’ in our plotting.

As rescaled_xda is still a Dask object, if we wanted to actually run the rescaling at this point (provided that all the data can fit into memory), we would use rescaled_xda.compute().

ML pipeline

Now that our data is in the proper shape and value range, we are ready to conduct spectral clustering. Here we will use a version of spectral clustering from dask_ml that is a scalable equivalent to operations from Scikit-learn that cluster pixels based on similarity (across all wavelength-bands, which makes it spectral clustering by spectra!)

client = Client(processes=False)
client

Client

Client-54b42401-2372-11ef-8aba-000d3ad44085

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: http://10.1.0.136:8787/status

Cluster Info

LocalCluster

ae410842

Dashboard: http://10.1.0.136:8787/status	Workers: 1
Total threads: 4	Total memory: 15.61 GiB
Status: running	Using processes: False

Scheduler Info

Scheduler

Scheduler-803cad1d-3321-43be-bde0-0c5c52a16326

Comm: inproc://10.1.0.136/2746/1	Workers: 1
Dashboard: http://10.1.0.136:8787/status	Total threads: 4
Started: Just now	Total memory: 15.61 GiB

Workers

Worker: 0

Comm: inproc://10.1.0.136/2746/4	Total threads: 4
Dashboard: http://10.1.0.136:35917/status	Memory: 15.61 GiB
Nanny: None
Local directory: /tmp/dask-scratch-space/worker-xkoxzjzy

Now we will compute and persist the rescaled data to feed into the ML pipeline. Notice that our X matrix below has the shape: n_samples, n_features as discussed earlier.

X = client.persist(rescaled_xda)
X.shape

(1071600, 8)

First we will set up the model with the number of clusters, and other options.

clf = SpectralClustering(
    n_clusters=4,
    random_state=0,
    gamma=None,
    kmeans_params={"init_max_iter": 5},
    persist_embedding=True,
)

This next step is the slow part. We’ll fit the model to our matrix X. Depending on your setup, it could take seconds to minutes to run depending on the size of our data.

%time clf.fit(X)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/distributed/client.py:3169: UserWarning: Sending large graph of size 81.80 MiB.
This may cause some slowdown.
Consider scattering data ahead of time and using futures.
  warnings.warn(

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/dask/base.py:1462: UserWarning: Running on a single-machine scheduler when a distributed client is active might lead to unexpected results.
  warnings.warn(

CPU times: user 23.8 s, sys: 2.86 s, total: 26.6 s
Wall time: 28.1 s

SpectralClustering(gamma=None, kmeans_params={'init_max_iter': 5}, n_clusters=4,
                   persist_embedding=True, random_state=0)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Let’s check the shape of the result:

labels = clf.assign_labels_.labels_.compute()
labels.shape

(1071600,)

labels

array([0, 2, 2, ..., 2, 2, 2], dtype=int32)

The result is a single vector of cluster labels.

Spectral Clustering for 1988

We have conducted the spectral clustering for 2017 and now we want to compare this result to the lake in 1988. Let’s load data from 1988 and run the same analysis as above.

We will use the same catalog, but we will search it for a different point in time and different Landsat mission

Load the data

bbox = [-118.89, 38.54, -118.57, 38.84]  # Region over a lake in Nevada, USA
datetime = "1988-06-01/1988-09-30"  # Summer months of 1988
collection = "landsat-c2-l2"
platform = "landsat-5"  # Searching through an earlier landsat mission
cloudy_less_than = 1  # percent

search = catalog.search(
    collections=["landsat-c2-l2"],
    bbox=bbox,
    datetime=datetime,
    query={"eo:cloud_cover": {"lt": cloudy_less_than}, "platform": {"in": [platform]}},
)

items = search.get_all_items()
item = items[1]  # select one of the results

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/pystac_client/item_search.py:849: FutureWarning: get_all_items() is deprecated, use item_collection() instead.
  warnings.warn(

Notice that Landsat 5 data from 1988 has slightly different spectra than Landsat 8 from 2017. Details like this are important to keep in mind when performing analyses that directly compare across missions.

assets = []
for _, asset in item.assets.items():
    try:
        assets.append(asset.extra_fields["eo:bands"][0])
    except:
        pass

cols_ordered = [
    "common_name",
    "description",
    "name",
    "center_wavelength",
    "full_width_half_max",
]
bands = pd.DataFrame.from_dict(assets)[cols_ordered]
bands

	common_name	description	name	center_wavelength	full_width_half_max
0	red	Visible red	TM_B3	0.66	0.06
1	blue	Visible blue	TM_B1	0.49	0.07
2	lwir	Long-wave infrared	TM_B6	11.45	2.10
3	green	Visible green	TM_B2	0.56	0.08
4	nir08	Near infrared	TM_B4	0.83	0.14
5	swir16	Short-wave infrared	TM_B5	1.65	0.20
6	swir22	Short-wave infrared	TM_B7	2.22	0.27

ds_1988 = odc.stac.stac_load(
    [item],
    bands=bands.common_name.values,
    bbox=bbox,
    chunks={},  # <-- use Dask
).isel(time=0)

epsg = item.properties["proj:epsg"]
ds_1988.attrs["crs"] = f"epsg:{epsg}"

da_1988 = ds_1988.to_array(dim="band")
da_1988

<xarray.DataArray (band: 7, y: 1128, x: 950)>
dask.array<stack, shape=(7, 1128, 950), dtype=uint16, chunksize=(1, 1128, 950), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
  * x            (x) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.638e+05
    spatial_ref  int32 32611
    time         datetime64[ns] 1988-07-02T18:03:59.010013
  * band         (band) object 'red' 'blue' 'lwir' ... 'nir08' 'swir16' 'swir22'
Attributes:
    crs:      epsg:32611

Reshape and Standardize

flattened_xda = da_1988.stack(z=("x", "y"))
flattened_t_xda = flattened_xda.transpose("z", "band")
with xr.set_options(keep_attrs=True):
    rescaled_xda = (flattened_t_xda - flattened_t_xda.mean()) / flattened_t_xda.std()
rescaled_xda

<xarray.DataArray (z: 1071600, band: 7)>
dask.array<truediv, shape=(1071600, 7), dtype=float64, chunksize=(1071600, 1), chunktype=numpy.ndarray>
Coordinates:
    spatial_ref  int32 32611
    time         datetime64[ns] 1988-07-02T18:03:59.010013
  * band         (band) object 'red' 'blue' 'lwir' ... 'nir08' 'swir16' 'swir22'
  * z            (z) object MultiIndex
  * x            (z) float64 3.353e+05 3.353e+05 ... 3.638e+05 3.638e+05
  * y            (z) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
Attributes:
    crs:      epsg:32611

Spectral Clustering

X = client.persist(rescaled_xda)
clf = SpectralClustering(
    n_clusters=4,
    random_state=0,
    gamma=None,
    kmeans_params={"init_max_iter": 5},
    persist_embedding=True,
)

%time clf.fit(X)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/distributed/client.py:3169: UserWarning: Sending large graph of size 73.62 MiB.
This may cause some slowdown.
Consider scattering data ahead of time and using futures.
  warnings.warn(

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/dask/base.py:1462: UserWarning: Running on a single-machine scheduler when a distributed client is active might lead to unexpected results.
  warnings.warn(

CPU times: user 19.5 s, sys: 2.53 s, total: 22.1 s
Wall time: 26.3 s

SpectralClustering(gamma=None, kmeans_params={'init_max_iter': 5}, n_clusters=4,
                   persist_embedding=True, random_state=0)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

labels = clf.assign_labels_.labels_.compute()
labels.shape

(1071600,)

labels

array([2, 0, 0, ..., 2, 2, 2], dtype=int32)

Spectral Clustering Over Time

Our hypothesis is that the lake’s area is receding over time and so we want to visualize the potential change. Let’s first visually compare the plot of the clustering results from the different time points.

result_plot_1988 + result_plot_2017

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/IPython/core/formatters.py:977, in MimeBundleFormatter.__call__(self, obj, include, exclude)
    974     method = get_real_method(obj, self.print_method)
    976     if method is not None:
--> 977         return method(include=include, exclude=exclude)
    978     return None
    979 else:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/core/dimension.py:1286, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1279 def _repr_mimebundle_(self, include=None, exclude=None):
   1280     """
   1281     Resolves the class hierarchy for the class rendering the
   1282     object using any display hooks registered on Store.display
   1283     hooks.  The output of all registered display_hooks is then
   1284     combined and returned.
   1285     """
-> 1286     return Store.render(self)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/core/options.py:1428, in Store.render(cls, obj)
   1426 data, metadata = {}, {}
   1427 for hook in hooks:
-> 1428     ret = hook(obj)
   1429     if ret is None:
   1430         continue

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:287, in pprint_display(obj)
    285 if not ip.display_formatter.formatters['text/plain'].pprint:
    286     return None
--> 287 return display(obj, raw_output=True)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:258, in display(obj, raw_output, **kwargs)
    256 elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    257     with option_state(obj):
--> 258         output = layout_display(obj)
    259 elif isinstance(obj, (HoloMap, DynamicMap)):
    260     with option_state(obj):

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:149, in display_hook.<locals>.wrapped(element)
    147 try:
    148     max_frames = OutputSettings.options['max_frames']
--> 149     mimebundle = fn(element, max_frames=max_frames)
    150     if mimebundle is None:
    151         return {}, {}

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:223, in layout_display(layout, max_frames)
    220     max_frame_warning(max_frames)
    221     return None
--> 223 return render(layout)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:76, in render(obj, **kwargs)
     73 if renderer.fig == 'pdf':
     74     renderer = renderer.instance(fig='png')
---> 76 return renderer.components(obj, **kwargs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:396, in Renderer.components(self, obj, fmt, comm, **kwargs)
    393 embed = (not (dynamic or streams or self.widget_mode == 'live') or config.embed)
    395 if embed or config.comms == 'default':
--> 396     return self._render_panel(plot, embed, comm)
    397 return self._render_ipywidget(plot)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:403, in Renderer._render_panel(self, plot, embed, comm)
    401 doc = Document()
    402 with config.set(embed=embed):
--> 403     model = plot.layout._render_model(doc, comm)
    404 if embed:
    405     return render_model(model, comm)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:748, in Viewable._render_model(self, doc, comm)
    746 if comm is None:
    747     comm = state._comm_manager.get_server_comm()
--> 748 model = self.get_root(doc, comm)
    750 if self._design and self._design.theme.bokeh_theme:
    751     doc.theme = self._design.theme.bokeh_theme

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:311, in Panel.get_root(self, doc, comm, preprocess)
    307 def get_root(
    308     self, doc: Optional[Document] = None, comm: Optional[Comm] = None,
    309     preprocess: bool = True
    310 ) -> Model:
--> 311     root = super().get_root(doc, comm, preprocess)
    312     # ALERT: Find a better way to handle this
    313     if hasattr(root, 'styles') and 'overflow-x' in root.styles:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:666, in Renderable.get_root(self, doc, comm, preprocess)
    664 wrapper = self._design._wrapper(self)
    665 if wrapper is self:
--> 666     root = self._get_model(doc, comm=comm)
    667     if preprocess:
    668         self._preprocess(root)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:177, in Panel._get_model(self, doc, root, parent, comm)
    175 root = root or model
    176 self._models[root.ref['id']] = (model, parent)
--> 177 objects, _ = self._get_objects(model, [], doc, root, comm)
    178 props = self._get_properties(doc)
    179 props[self._property_mapping['objects']] = objects

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:159, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    157 else:
    158     try:
--> 159         child = pane._get_model(doc, root, model, comm)
    160     except RerenderError as e:
    161         if e.layout is not None and e.layout is not self:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:423, in HoloViews._get_model(self, doc, root, parent, comm)
    421     plot = self.object
    422 else:
--> 423     plot = self._render(doc, comm, root)
    425 plot.pane = self
    426 backend = plot.renderer.backend

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:518, in HoloViews._render(self, doc, comm, root)
    515     if comm:
    516         kwargs['comm'] = comm
--> 518 return renderer.get_plot(self.object, **kwargs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/renderer.py:68, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     61 @bothmethod
     62 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     63     """
     64     Given a HoloViews Viewable return a corresponding plot instance.
     65     Allows supplying a document attach the plot to, useful when
     66     combining the bokeh model with another plot.
     67     """
---> 68     plot = super().get_plot(obj, doc, renderer, **kwargs)
     69     if plot.document is None:
     70         plot.document = Document() if self_or_cls.notebook_context else curdoc()

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:240, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    237     defaults = [kd.default for kd in plot.dimensions]
    238     init_key = tuple(v if d is None else d for v, d in
    239                      zip(plot.keys[0], defaults))
--> 240     plot.update(init_key)
    241 else:
    242     plot = obj

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/plot.py:955, in DimensionedPlot.update(self, key)
    953 def update(self, key):
    954     if len(self) == 1 and key in (0, self.keys[0]) and not self.drawn:
--> 955         return self.initialize_plot()
    956     item = self.__getitem__(key)
    957     self.traverse(lambda x: setattr(x, '_updated', True))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:945, in LayoutPlot.initialize_plot(self, plots, ranges)
    942     continue
    944 shared_plots = list(passed_plots) if self.shared_axes else None
--> 945 subplots = subplot.initialize_plot(ranges=ranges, plots=shared_plots)
    946 nsubplots = len(subplots)
    948 modes = {sp.sizing_mode for sp in subplots
    949          if sp.sizing_mode not in (None, 'auto', 'fixed')}

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:1094, in AdjointLayoutPlot.initialize_plot(self, ranges, plots)
   1092     else:
   1093         passed_plots = plots + adjoined_plots
-> 1094         adjoined_plots.append(subplot.initialize_plot(ranges=ranges, plots=passed_plots))
   1095 self.drawn = True
   1096 if not adjoined_plots: adjoined_plots = [None]

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/bokeh/plot.py:108, in GeoPlot.initialize_plot(self, ranges, plot, plots, source)
    106 def initialize_plot(self, ranges=None, plot=None, plots=None, source=None):
    107     opts = {} if isinstance(self, HvOverlayPlot) else {'source': source}
--> 108     fig = super().initialize_plot(ranges, plot, plots, **opts)
    109     if self.geographic and self.show_bounds and not self.overlaid:
    110         from . import GeoShapePlot

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:1875, in ElementPlot.initialize_plot(self, ranges, plot, plots, source)
   1873 # Initialize plot, source and glyph
   1874 if plot is None:
-> 1875     plot = self._init_plot(key, style_element, ranges=ranges, plots=plots)
   1876     self._populate_axis_handles(plot)
   1877 else:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:644, in ElementPlot._init_plot(self, key, element, plots, ranges)
    641 subplots = list(self.subplots.values()) if self.subplots else []
    643 axis_specs = {'x': {}, 'y': {}}
--> 644 axis_specs['x']['x'] = self._axis_props(plots, subplots, element, ranges, pos=0) + (self.xaxis, {})
    645 if self.multi_y:
    646     if not bokeh32:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:477, in ElementPlot._axis_props(self, plots, subplots, element, ranges, pos, dim, range_tags_extras, extra_range_name)
    475 else:
    476     try:
--> 477         l, b, r, t = self.get_extents(range_el, ranges, dimension=dim)
    478     except TypeError:
    479         # Backward compatibility for e.g. GeoViews=<1.10.1 since dimension
    480         # is a newly added keyword argument in HoloViews 1.17
    481         l, b, r, t = self.get_extents(range_el, ranges)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/plot.py:73, in ProjectionPlot.get_extents(self, element, ranges, range_type, **kwargs)
     71     extents = None
     72 else:
---> 73     extents = project_extents(extents, element.crs, proj)
     74 return (np.nan,)*4 if not extents else extents

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/util.py:100, in project_extents(extents, src_proj, dest_proj, tol)
     98     geom_in_src_proj = geom_clipped_to_dest_proj
     99 try:
--> 100     geom_in_crs = dest_proj.project_geometry(geom_in_src_proj, src_proj)
    101 except ValueError as e:
    102     src_name =type(src_proj).__name__

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:824, in Projection.project_geometry(self, geometry, src_crs)
    822 if not method_name:
    823     raise ValueError(f'Unsupported geometry type {geom_type!r}')
--> 824 return getattr(self, method_name)(geometry, src_crs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:960, in Projection._project_polygon(self, polygon, src_crs)
    958     is_ccw = True
    959 else:
--> 960     is_ccw = polygon.exterior.is_ccw
    961 # Project the polygon exterior/interior rings.
    962 # Each source ring will result in either a ring, or one or more
    963 # lines.
    964 rings = []

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/geometry/polygon.py:99, in LinearRing.is_ccw(self)
     96 @property
     97 def is_ccw(self):
     98     """True is the ring is oriented counter clock-wise"""
---> 99     return bool(self.impl['is_ccw'](self))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/algorithms/cga.py:14, in is_ccw_impl.<locals>.is_ccw_op(ring)
     13 def is_ccw_op(ring):
---> 14     return signed_area(ring) >= 0.0

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/algorithms/cga.py:7, in signed_area(ring)
      3 """Return the signed area enclosed by a ring in linear time using the
      4 algorithm at: https://web.archive.org/web/20080209143651/http://cgafaq.info:80/wiki/Polygon_Area
      5 """
      6 xs, ys = ring.coords.xy
----> 7 xs.append(xs[1])
      8 ys.append(ys[1])
      9 return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(ring.coords)))/2.0

IndexError: array index out of range

:Layout
   .Image.I  :Image   [x,y]   (value)
   .Image.II :Image   [x,y]   (value)

By hovering over the lake in each image, we can see that the water was labeled (‘value’) with a certain cluster number in both images. We will programmatically grab the water cluster value from the middle of the lake using pyproj to convert from longtitude/latitude coordinates.

lon_lake_center = -118.71
lat_lake_center = 38.7

proj = pyproj.Proj(unstacked_1988.crs)
lake_center_x, lake_center_y = proj(lon_lake_center, lat_lake_center)

water_cluster_1988 = int(unstacked_1988.sel(x=lake_center_x, y=lake_center_y, method='nearest'))
water_cluster_2017 = int(unstacked_2017.sel(x=lake_center_x, y=lake_center_y, method='nearest'))

print('water cluster values:', water_cluster_1988, water_cluster_2017)

water cluster values: 3 1

Now, let’s set any value that isn’t our water cluster label to 0.

with xr.set_options(keep_attrs=True):
    water_1988 = (unstacked_1988 == water_cluster_1988).astype(int)
    water_2017 = (unstacked_2017 == water_cluster_2017).astype(int)

water_1988_plot = water_1988.hvplot(
    x="x", y="y", cmap="greys", geo=True, colorbar=False, title="1988 Water"
)

water_2017_plot = water_2017.hvplot(
    x="x", y="y", cmap="greys", geo=True, colorbar=False, title="2017 Water"
)

water_1988_plot + water_2017_plot

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/IPython/core/formatters.py:977, in MimeBundleFormatter.__call__(self, obj, include, exclude)
    974     method = get_real_method(obj, self.print_method)
    976     if method is not None:
--> 977         return method(include=include, exclude=exclude)
    978     return None
    979 else:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/core/dimension.py:1286, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1279 def _repr_mimebundle_(self, include=None, exclude=None):
   1280     """
   1281     Resolves the class hierarchy for the class rendering the
   1282     object using any display hooks registered on Store.display
   1283     hooks.  The output of all registered display_hooks is then
   1284     combined and returned.
   1285     """
-> 1286     return Store.render(self)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/core/options.py:1428, in Store.render(cls, obj)
   1426 data, metadata = {}, {}
   1427 for hook in hooks:
-> 1428     ret = hook(obj)
   1429     if ret is None:
   1430         continue

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:287, in pprint_display(obj)
    285 if not ip.display_formatter.formatters['text/plain'].pprint:
    286     return None
--> 287 return display(obj, raw_output=True)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:258, in display(obj, raw_output, **kwargs)
    256 elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    257     with option_state(obj):
--> 258         output = layout_display(obj)
    259 elif isinstance(obj, (HoloMap, DynamicMap)):
    260     with option_state(obj):

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:149, in display_hook.<locals>.wrapped(element)
    147 try:
    148     max_frames = OutputSettings.options['max_frames']
--> 149     mimebundle = fn(element, max_frames=max_frames)
    150     if mimebundle is None:
    151         return {}, {}

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:223, in layout_display(layout, max_frames)
    220     max_frame_warning(max_frames)
    221     return None
--> 223 return render(layout)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:76, in render(obj, **kwargs)
     73 if renderer.fig == 'pdf':
     74     renderer = renderer.instance(fig='png')
---> 76 return renderer.components(obj, **kwargs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:396, in Renderer.components(self, obj, fmt, comm, **kwargs)
    393 embed = (not (dynamic or streams or self.widget_mode == 'live') or config.embed)
    395 if embed or config.comms == 'default':
--> 396     return self._render_panel(plot, embed, comm)
    397 return self._render_ipywidget(plot)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:403, in Renderer._render_panel(self, plot, embed, comm)
    401 doc = Document()
    402 with config.set(embed=embed):
--> 403     model = plot.layout._render_model(doc, comm)
    404 if embed:
    405     return render_model(model, comm)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:748, in Viewable._render_model(self, doc, comm)
    746 if comm is None:
    747     comm = state._comm_manager.get_server_comm()
--> 748 model = self.get_root(doc, comm)
    750 if self._design and self._design.theme.bokeh_theme:
    751     doc.theme = self._design.theme.bokeh_theme

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:311, in Panel.get_root(self, doc, comm, preprocess)
    307 def get_root(
    308     self, doc: Optional[Document] = None, comm: Optional[Comm] = None,
    309     preprocess: bool = True
    310 ) -> Model:
--> 311     root = super().get_root(doc, comm, preprocess)
    312     # ALERT: Find a better way to handle this
    313     if hasattr(root, 'styles') and 'overflow-x' in root.styles:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:666, in Renderable.get_root(self, doc, comm, preprocess)
    664 wrapper = self._design._wrapper(self)
    665 if wrapper is self:
--> 666     root = self._get_model(doc, comm=comm)
    667     if preprocess:
    668         self._preprocess(root)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:177, in Panel._get_model(self, doc, root, parent, comm)
    175 root = root or model
    176 self._models[root.ref['id']] = (model, parent)
--> 177 objects, _ = self._get_objects(model, [], doc, root, comm)
    178 props = self._get_properties(doc)
    179 props[self._property_mapping['objects']] = objects

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:159, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    157 else:
    158     try:
--> 159         child = pane._get_model(doc, root, model, comm)
    160     except RerenderError as e:
    161         if e.layout is not None and e.layout is not self:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:423, in HoloViews._get_model(self, doc, root, parent, comm)
    421     plot = self.object
    422 else:
--> 423     plot = self._render(doc, comm, root)
    425 plot.pane = self
    426 backend = plot.renderer.backend

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:518, in HoloViews._render(self, doc, comm, root)
    515     if comm:
    516         kwargs['comm'] = comm
--> 518 return renderer.get_plot(self.object, **kwargs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/renderer.py:68, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     61 @bothmethod
     62 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     63     """
     64     Given a HoloViews Viewable return a corresponding plot instance.
     65     Allows supplying a document attach the plot to, useful when
     66     combining the bokeh model with another plot.
     67     """
---> 68     plot = super().get_plot(obj, doc, renderer, **kwargs)
     69     if plot.document is None:
     70         plot.document = Document() if self_or_cls.notebook_context else curdoc()

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:240, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    237     defaults = [kd.default for kd in plot.dimensions]
    238     init_key = tuple(v if d is None else d for v, d in
    239                      zip(plot.keys[0], defaults))
--> 240     plot.update(init_key)
    241 else:
    242     plot = obj

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/plot.py:955, in DimensionedPlot.update(self, key)
    953 def update(self, key):
    954     if len(self) == 1 and key in (0, self.keys[0]) and not self.drawn:
--> 955         return self.initialize_plot()
    956     item = self.__getitem__(key)
    957     self.traverse(lambda x: setattr(x, '_updated', True))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:945, in LayoutPlot.initialize_plot(self, plots, ranges)
    942     continue
    944 shared_plots = list(passed_plots) if self.shared_axes else None
--> 945 subplots = subplot.initialize_plot(ranges=ranges, plots=shared_plots)
    946 nsubplots = len(subplots)
    948 modes = {sp.sizing_mode for sp in subplots
    949          if sp.sizing_mode not in (None, 'auto', 'fixed')}

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:1094, in AdjointLayoutPlot.initialize_plot(self, ranges, plots)
   1092     else:
   1093         passed_plots = plots + adjoined_plots
-> 1094         adjoined_plots.append(subplot.initialize_plot(ranges=ranges, plots=passed_plots))
   1095 self.drawn = True
   1096 if not adjoined_plots: adjoined_plots = [None]

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/bokeh/plot.py:108, in GeoPlot.initialize_plot(self, ranges, plot, plots, source)
    106 def initialize_plot(self, ranges=None, plot=None, plots=None, source=None):
    107     opts = {} if isinstance(self, HvOverlayPlot) else {'source': source}
--> 108     fig = super().initialize_plot(ranges, plot, plots, **opts)
    109     if self.geographic and self.show_bounds and not self.overlaid:
    110         from . import GeoShapePlot

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:1875, in ElementPlot.initialize_plot(self, ranges, plot, plots, source)
   1873 # Initialize plot, source and glyph
   1874 if plot is None:
-> 1875     plot = self._init_plot(key, style_element, ranges=ranges, plots=plots)
   1876     self._populate_axis_handles(plot)
   1877 else:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:644, in ElementPlot._init_plot(self, key, element, plots, ranges)
    641 subplots = list(self.subplots.values()) if self.subplots else []
    643 axis_specs = {'x': {}, 'y': {}}
--> 644 axis_specs['x']['x'] = self._axis_props(plots, subplots, element, ranges, pos=0) + (self.xaxis, {})
    645 if self.multi_y:
    646     if not bokeh32:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:477, in ElementPlot._axis_props(self, plots, subplots, element, ranges, pos, dim, range_tags_extras, extra_range_name)
    475 else:
    476     try:
--> 477         l, b, r, t = self.get_extents(range_el, ranges, dimension=dim)
    478     except TypeError:
    479         # Backward compatibility for e.g. GeoViews=<1.10.1 since dimension
    480         # is a newly added keyword argument in HoloViews 1.17
    481         l, b, r, t = self.get_extents(range_el, ranges)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/plot.py:73, in ProjectionPlot.get_extents(self, element, ranges, range_type, **kwargs)
     71     extents = None
     72 else:
---> 73     extents = project_extents(extents, element.crs, proj)
     74 return (np.nan,)*4 if not extents else extents

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/util.py:100, in project_extents(extents, src_proj, dest_proj, tol)
     98     geom_in_src_proj = geom_clipped_to_dest_proj
     99 try:
--> 100     geom_in_crs = dest_proj.project_geometry(geom_in_src_proj, src_proj)
    101 except ValueError as e:
    102     src_name =type(src_proj).__name__

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:824, in Projection.project_geometry(self, geometry, src_crs)
    822 if not method_name:
    823     raise ValueError(f'Unsupported geometry type {geom_type!r}')
--> 824 return getattr(self, method_name)(geometry, src_crs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:960, in Projection._project_polygon(self, polygon, src_crs)
    958     is_ccw = True
    959 else:
--> 960     is_ccw = polygon.exterior.is_ccw
    961 # Project the polygon exterior/interior rings.
    962 # Each source ring will result in either a ring, or one or more
    963 # lines.
    964 rings = []

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/geometry/polygon.py:99, in LinearRing.is_ccw(self)
     96 @property
     97 def is_ccw(self):
     98     """True is the ring is oriented counter clock-wise"""
---> 99     return bool(self.impl['is_ccw'](self))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/algorithms/cga.py:14, in is_ccw_impl.<locals>.is_ccw_op(ring)
     13 def is_ccw_op(ring):
---> 14     return signed_area(ring) >= 0.0

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/algorithms/cga.py:7, in signed_area(ring)
      3 """Return the signed area enclosed by a ring in linear time using the
      4 algorithm at: https://web.archive.org/web/20080209143651/http://cgafaq.info:80/wiki/Polygon_Area
      5 """
      6 xs, ys = ring.coords.xy
----> 7 xs.append(xs[1])
      8 ys.append(ys[1])
      9 return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(ring.coords)))/2.0

IndexError: array index out of range

:Layout
   .Image.I  :Image   [x,y]   (value)
   .Image.II :Image   [x,y]   (value)

Now we can take the difference of these water label arrays to see exactly where the water levels has changed.

with xr.set_options(keep_attrs=True):
    water_diff = water_1988 - water_2017

Red pixels (array value ‘1’) of our image below are where water was lost from 1988 to 2017.

water_diff.hvplot(
    x="x", y="y", cmap='coolwarm', geo=True, xlabel="long", ylabel="lat", colorbar=False, title="Water Change 1988-2017",
)

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/IPython/core/formatters.py:977, in MimeBundleFormatter.__call__(self, obj, include, exclude)
    974     method = get_real_method(obj, self.print_method)
    976     if method is not None:
--> 977         return method(include=include, exclude=exclude)
    978     return None
    979 else:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/core/dimension.py:1286, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1279 def _repr_mimebundle_(self, include=None, exclude=None):
   1280     """
   1281     Resolves the class hierarchy for the class rendering the
   1282     object using any display hooks registered on Store.display
   1283     hooks.  The output of all registered display_hooks is then
   1284     combined and returned.
   1285     """
-> 1286     return Store.render(self)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/core/options.py:1428, in Store.render(cls, obj)
   1426 data, metadata = {}, {}
   1427 for hook in hooks:
-> 1428     ret = hook(obj)
   1429     if ret is None:
   1430         continue

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:287, in pprint_display(obj)
    285 if not ip.display_formatter.formatters['text/plain'].pprint:
    286     return None
--> 287 return display(obj, raw_output=True)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:255, in display(obj, raw_output, **kwargs)
    253 elif isinstance(obj, (CompositeOverlay, ViewableElement)):
    254     with option_state(obj):
--> 255         output = element_display(obj)
    256 elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    257     with option_state(obj):

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:149, in display_hook.<locals>.wrapped(element)
    147 try:
    148     max_frames = OutputSettings.options['max_frames']
--> 149     mimebundle = fn(element, max_frames=max_frames)
    150     if mimebundle is None:
    151         return {}, {}

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:195, in element_display(element, max_frames)
    192 if type(element) not in Store.registry[backend]:
    193     return None
--> 195 return render(element)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:76, in render(obj, **kwargs)
     73 if renderer.fig == 'pdf':
     74     renderer = renderer.instance(fig='png')
---> 76 return renderer.components(obj, **kwargs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:396, in Renderer.components(self, obj, fmt, comm, **kwargs)
    393 embed = (not (dynamic or streams or self.widget_mode == 'live') or config.embed)
    395 if embed or config.comms == 'default':
--> 396     return self._render_panel(plot, embed, comm)
    397 return self._render_ipywidget(plot)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:403, in Renderer._render_panel(self, plot, embed, comm)
    401 doc = Document()
    402 with config.set(embed=embed):
--> 403     model = plot.layout._render_model(doc, comm)
    404 if embed:
    405     return render_model(model, comm)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:748, in Viewable._render_model(self, doc, comm)
    746 if comm is None:
    747     comm = state._comm_manager.get_server_comm()
--> 748 model = self.get_root(doc, comm)
    750 if self._design and self._design.theme.bokeh_theme:
    751     doc.theme = self._design.theme.bokeh_theme

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:311, in Panel.get_root(self, doc, comm, preprocess)
    307 def get_root(
    308     self, doc: Optional[Document] = None, comm: Optional[Comm] = None,
    309     preprocess: bool = True
    310 ) -> Model:
--> 311     root = super().get_root(doc, comm, preprocess)
    312     # ALERT: Find a better way to handle this
    313     if hasattr(root, 'styles') and 'overflow-x' in root.styles:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:666, in Renderable.get_root(self, doc, comm, preprocess)
    664 wrapper = self._design._wrapper(self)
    665 if wrapper is self:
--> 666     root = self._get_model(doc, comm=comm)
    667     if preprocess:
    668         self._preprocess(root)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:177, in Panel._get_model(self, doc, root, parent, comm)
    175 root = root or model
    176 self._models[root.ref['id']] = (model, parent)
--> 177 objects, _ = self._get_objects(model, [], doc, root, comm)
    178 props = self._get_properties(doc)
    179 props[self._property_mapping['objects']] = objects

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:159, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    157 else:
    158     try:
--> 159         child = pane._get_model(doc, root, model, comm)
    160     except RerenderError as e:
    161         if e.layout is not None and e.layout is not self:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:423, in HoloViews._get_model(self, doc, root, parent, comm)
    421     plot = self.object
    422 else:
--> 423     plot = self._render(doc, comm, root)
    425 plot.pane = self
    426 backend = plot.renderer.backend

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:518, in HoloViews._render(self, doc, comm, root)
    515     if comm:
    516         kwargs['comm'] = comm
--> 518 return renderer.get_plot(self.object, **kwargs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/renderer.py:68, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     61 @bothmethod
     62 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     63     """
     64     Given a HoloViews Viewable return a corresponding plot instance.
     65     Allows supplying a document attach the plot to, useful when
     66     combining the bokeh model with another plot.
     67     """
---> 68     plot = super().get_plot(obj, doc, renderer, **kwargs)
     69     if plot.document is None:
     70         plot.document = Document() if self_or_cls.notebook_context else curdoc()

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:240, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    237     defaults = [kd.default for kd in plot.dimensions]
    238     init_key = tuple(v if d is None else d for v, d in
    239                      zip(plot.keys[0], defaults))
--> 240     plot.update(init_key)
    241 else:
    242     plot = obj

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/plot.py:955, in DimensionedPlot.update(self, key)
    953 def update(self, key):
    954     if len(self) == 1 and key in (0, self.keys[0]) and not self.drawn:
--> 955         return self.initialize_plot()
    956     item = self.__getitem__(key)
    957     self.traverse(lambda x: setattr(x, '_updated', True))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/bokeh/plot.py:108, in GeoPlot.initialize_plot(self, ranges, plot, plots, source)
    106 def initialize_plot(self, ranges=None, plot=None, plots=None, source=None):
    107     opts = {} if isinstance(self, HvOverlayPlot) else {'source': source}
--> 108     fig = super().initialize_plot(ranges, plot, plots, **opts)
    109     if self.geographic and self.show_bounds and not self.overlaid:
    110         from . import GeoShapePlot

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:1875, in ElementPlot.initialize_plot(self, ranges, plot, plots, source)
   1873 # Initialize plot, source and glyph
   1874 if plot is None:
-> 1875     plot = self._init_plot(key, style_element, ranges=ranges, plots=plots)
   1876     self._populate_axis_handles(plot)
   1877 else:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:644, in ElementPlot._init_plot(self, key, element, plots, ranges)
    641 subplots = list(self.subplots.values()) if self.subplots else []
    643 axis_specs = {'x': {}, 'y': {}}
--> 644 axis_specs['x']['x'] = self._axis_props(plots, subplots, element, ranges, pos=0) + (self.xaxis, {})
    645 if self.multi_y:
    646     if not bokeh32:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:477, in ElementPlot._axis_props(self, plots, subplots, element, ranges, pos, dim, range_tags_extras, extra_range_name)
    475 else:
    476     try:
--> 477         l, b, r, t = self.get_extents(range_el, ranges, dimension=dim)
    478     except TypeError:
    479         # Backward compatibility for e.g. GeoViews=<1.10.1 since dimension
    480         # is a newly added keyword argument in HoloViews 1.17
    481         l, b, r, t = self.get_extents(range_el, ranges)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/plot.py:73, in ProjectionPlot.get_extents(self, element, ranges, range_type, **kwargs)
     71     extents = None
     72 else:
---> 73     extents = project_extents(extents, element.crs, proj)
     74 return (np.nan,)*4 if not extents else extents

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/util.py:100, in project_extents(extents, src_proj, dest_proj, tol)
     98     geom_in_src_proj = geom_clipped_to_dest_proj
     99 try:
--> 100     geom_in_crs = dest_proj.project_geometry(geom_in_src_proj, src_proj)
    101 except ValueError as e:
    102     src_name =type(src_proj).__name__

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:824, in Projection.project_geometry(self, geometry, src_crs)
    822 if not method_name:
    823     raise ValueError(f'Unsupported geometry type {geom_type!r}')
--> 824 return getattr(self, method_name)(geometry, src_crs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:960, in Projection._project_polygon(self, polygon, src_crs)
    958     is_ccw = True
    959 else:
--> 960     is_ccw = polygon.exterior.is_ccw
    961 # Project the polygon exterior/interior rings.
    962 # Each source ring will result in either a ring, or one or more
    963 # lines.
    964 rings = []

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/geometry/polygon.py:99, in LinearRing.is_ccw(self)
     96 @property
     97 def is_ccw(self):
     98     """True is the ring is oriented counter clock-wise"""
---> 99     return bool(self.impl['is_ccw'](self))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/algorithms/cga.py:14, in is_ccw_impl.<locals>.is_ccw_op(ring)
     13 def is_ccw_op(ring):
---> 14     return signed_area(ring) >= 0.0

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/algorithms/cga.py:7, in signed_area(ring)
      3 """Return the signed area enclosed by a ring in linear time using the
      4 algorithm at: https://web.archive.org/web/20080209143651/http://cgafaq.info:80/wiki/Polygon_Area
      5 """
      6 xs, ys = ring.coords.xy
----> 7 xs.append(xs[1])
      8 ys.append(ys[1])
      9 return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(ring.coords)))/2.0

IndexError: array index out of range

:Image   [x,y]   (value)

We did it! We are observing the change in the lake shoreline over time using a simple spectral clustering approach.

Let’s finish things off by adding some geo tiles as a background. To only display the colored pixels overlaid on geo tiles, we could either set the array’s background value (‘0’) to ‘Not a Number’ (NaN), or we could just inform hvPlot that we want the background valued pixels to be transparent with .redim.nodata(value=0).

water_diff.hvplot(
        x="x", y="y", width=400, height=400, cmap='coolwarm', geo=True, xlabel="lon", ylabel="lat", alpha=1, colorbar=False, title="Water Loss from 1988 to 2017", tiles="ESRI",
).redim.nodata(value=0)

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/IPython/core/formatters.py:977, in MimeBundleFormatter.__call__(self, obj, include, exclude)
    974     method = get_real_method(obj, self.print_method)
    976     if method is not None:
--> 977         return method(include=include, exclude=exclude)
    978     return None
    979 else:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/core/dimension.py:1286, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1279 def _repr_mimebundle_(self, include=None, exclude=None):
   1280     """
   1281     Resolves the class hierarchy for the class rendering the
   1282     object using any display hooks registered on Store.display
   1283     hooks.  The output of all registered display_hooks is then
   1284     combined and returned.
   1285     """
-> 1286     return Store.render(self)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/core/options.py:1428, in Store.render(cls, obj)
   1426 data, metadata = {}, {}
   1427 for hook in hooks:
-> 1428     ret = hook(obj)
   1429     if ret is None:
   1430         continue

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:287, in pprint_display(obj)
    285 if not ip.display_formatter.formatters['text/plain'].pprint:
    286     return None
--> 287 return display(obj, raw_output=True)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:255, in display(obj, raw_output, **kwargs)
    253 elif isinstance(obj, (CompositeOverlay, ViewableElement)):
    254     with option_state(obj):
--> 255         output = element_display(obj)
    256 elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    257     with option_state(obj):

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:149, in display_hook.<locals>.wrapped(element)
    147 try:
    148     max_frames = OutputSettings.options['max_frames']
--> 149     mimebundle = fn(element, max_frames=max_frames)
    150     if mimebundle is None:
    151         return {}, {}

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:195, in element_display(element, max_frames)
    192 if type(element) not in Store.registry[backend]:
    193     return None
--> 195 return render(element)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/ipython/display_hooks.py:76, in render(obj, **kwargs)
     73 if renderer.fig == 'pdf':
     74     renderer = renderer.instance(fig='png')
---> 76 return renderer.components(obj, **kwargs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:396, in Renderer.components(self, obj, fmt, comm, **kwargs)
    393 embed = (not (dynamic or streams or self.widget_mode == 'live') or config.embed)
    395 if embed or config.comms == 'default':
--> 396     return self._render_panel(plot, embed, comm)
    397 return self._render_ipywidget(plot)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:403, in Renderer._render_panel(self, plot, embed, comm)
    401 doc = Document()
    402 with config.set(embed=embed):
--> 403     model = plot.layout._render_model(doc, comm)
    404 if embed:
    405     return render_model(model, comm)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:748, in Viewable._render_model(self, doc, comm)
    746 if comm is None:
    747     comm = state._comm_manager.get_server_comm()
--> 748 model = self.get_root(doc, comm)
    750 if self._design and self._design.theme.bokeh_theme:
    751     doc.theme = self._design.theme.bokeh_theme

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:311, in Panel.get_root(self, doc, comm, preprocess)
    307 def get_root(
    308     self, doc: Optional[Document] = None, comm: Optional[Comm] = None,
    309     preprocess: bool = True
    310 ) -> Model:
--> 311     root = super().get_root(doc, comm, preprocess)
    312     # ALERT: Find a better way to handle this
    313     if hasattr(root, 'styles') and 'overflow-x' in root.styles:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:666, in Renderable.get_root(self, doc, comm, preprocess)
    664 wrapper = self._design._wrapper(self)
    665 if wrapper is self:
--> 666     root = self._get_model(doc, comm=comm)
    667     if preprocess:
    668         self._preprocess(root)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:177, in Panel._get_model(self, doc, root, parent, comm)
    175 root = root or model
    176 self._models[root.ref['id']] = (model, parent)
--> 177 objects, _ = self._get_objects(model, [], doc, root, comm)
    178 props = self._get_properties(doc)
    179 props[self._property_mapping['objects']] = objects

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:159, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    157 else:
    158     try:
--> 159         child = pane._get_model(doc, root, model, comm)
    160     except RerenderError as e:
    161         if e.layout is not None and e.layout is not self:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:423, in HoloViews._get_model(self, doc, root, parent, comm)
    421     plot = self.object
    422 else:
--> 423     plot = self._render(doc, comm, root)
    425 plot.pane = self
    426 backend = plot.renderer.backend

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:518, in HoloViews._render(self, doc, comm, root)
    515     if comm:
    516         kwargs['comm'] = comm
--> 518 return renderer.get_plot(self.object, **kwargs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/renderer.py:68, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     61 @bothmethod
     62 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     63     """
     64     Given a HoloViews Viewable return a corresponding plot instance.
     65     Allows supplying a document attach the plot to, useful when
     66     combining the bokeh model with another plot.
     67     """
---> 68     plot = super().get_plot(obj, doc, renderer, **kwargs)
     69     if plot.document is None:
     70         plot.document = Document() if self_or_cls.notebook_context else curdoc()

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/renderer.py:240, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    237     defaults = [kd.default for kd in plot.dimensions]
    238     init_key = tuple(v if d is None else d for v, d in
    239                      zip(plot.keys[0], defaults))
--> 240     plot.update(init_key)
    241 else:
    242     plot = obj

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/plot.py:955, in DimensionedPlot.update(self, key)
    953 def update(self, key):
    954     if len(self) == 1 and key in (0, self.keys[0]) and not self.drawn:
--> 955         return self.initialize_plot()
    956     item = self.__getitem__(key)
    957     self.traverse(lambda x: setattr(x, '_updated', True))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/bokeh/plot.py:108, in GeoPlot.initialize_plot(self, ranges, plot, plots, source)
    106 def initialize_plot(self, ranges=None, plot=None, plots=None, source=None):
    107     opts = {} if isinstance(self, HvOverlayPlot) else {'source': source}
--> 108     fig = super().initialize_plot(ranges, plot, plots, **opts)
    109     if self.geographic and self.show_bounds and not self.overlaid:
    110         from . import GeoShapePlot

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:2886, in OverlayPlot.initialize_plot(self, ranges, plot, plots)
   2884 self.tabs = self.tabs or any(isinstance(sp, TablePlot) for sp in self.subplots.values())
   2885 if plot is None and not self.tabs and not self.batched:
-> 2886     plot = self._init_plot(key, element, ranges=ranges, plots=plots)
   2887     self._populate_axis_handles(plot)
   2888 self.handles['plot'] = plot

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:644, in ElementPlot._init_plot(self, key, element, plots, ranges)
    641 subplots = list(self.subplots.values()) if self.subplots else []
    643 axis_specs = {'x': {}, 'y': {}}
--> 644 axis_specs['x']['x'] = self._axis_props(plots, subplots, element, ranges, pos=0) + (self.xaxis, {})
    645 if self.multi_y:
    646     if not bokeh32:

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:477, in ElementPlot._axis_props(self, plots, subplots, element, ranges, pos, dim, range_tags_extras, extra_range_name)
    475 else:
    476     try:
--> 477         l, b, r, t = self.get_extents(range_el, ranges, dimension=dim)
    478     except TypeError:
    479         # Backward compatibility for e.g. GeoViews=<1.10.1 since dimension
    480         # is a newly added keyword argument in HoloViews 1.17
    481         l, b, r, t = self.get_extents(range_el, ranges)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/plot.py:67, in ProjectionPlot.get_extents(self, element, ranges, range_type, **kwargs)
     65     (x0, x1), (y0, y1) = proj.x_limits, proj.y_limits
     66     return (x0, y0, x1, y1)
---> 67 extents = super().get_extents(element, ranges, range_type)
     68 if not getattr(element, 'crs', None) or not self.geographic:
     69     return extents

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/plot.py:1986, in GenericOverlayPlot.get_extents(self, overlay, ranges, range_type, dimension, **kwargs)
   1985 def get_extents(self, overlay, ranges, range_type='combined', dimension=None, **kwargs):
-> 1986     subplot_extents = self._get_subplot_extents(overlay, ranges, range_type, dimension)
   1987     zrange = isinstance(self.projection, str) and self.projection == '3d'
   1988     extents = {k: util.max_extents(rs, zrange) for k, rs in subplot_extents.items()}

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/plot.py:1981, in GenericOverlayPlot._get_subplot_extents(self, overlay, ranges, range_type, dimension)
   1979         sp_ranges = util.match_spec(layer, ranges) if ranges else {}
   1980     for rt in extents:
-> 1981         extent = subplot.get_extents(layer, sp_ranges, range_type=rt)
   1982         extents[rt].append(extent)
   1983 return extents

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/plot.py:73, in ProjectionPlot.get_extents(self, element, ranges, range_type, **kwargs)
     71     extents = None
     72 else:
---> 73     extents = project_extents(extents, element.crs, proj)
     74 return (np.nan,)*4 if not extents else extents

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/util.py:100, in project_extents(extents, src_proj, dest_proj, tol)
     98     geom_in_src_proj = geom_clipped_to_dest_proj
     99 try:
--> 100     geom_in_crs = dest_proj.project_geometry(geom_in_src_proj, src_proj)
    101 except ValueError as e:
    102     src_name =type(src_proj).__name__

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:824, in Projection.project_geometry(self, geometry, src_crs)
    822 if not method_name:
    823     raise ValueError(f'Unsupported geometry type {geom_type!r}')
--> 824 return getattr(self, method_name)(geometry, src_crs)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:960, in Projection._project_polygon(self, polygon, src_crs)
    958     is_ccw = True
    959 else:
--> 960     is_ccw = polygon.exterior.is_ccw
    961 # Project the polygon exterior/interior rings.
    962 # Each source ring will result in either a ring, or one or more
    963 # lines.
    964 rings = []

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/geometry/polygon.py:99, in LinearRing.is_ccw(self)
     96 @property
     97 def is_ccw(self):
     98     """True is the ring is oriented counter clock-wise"""
---> 99     return bool(self.impl['is_ccw'](self))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/algorithms/cga.py:14, in is_ccw_impl.<locals>.is_ccw_op(ring)
     13 def is_ccw_op(ring):
---> 14     return signed_area(ring) >= 0.0

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/shapely/algorithms/cga.py:7, in signed_area(ring)
      3 """Return the signed area enclosed by a ring in linear time using the
      4 algorithm at: https://web.archive.org/web/20080209143651/http://cgafaq.info:80/wiki/Polygon_Area
      5 """
      6 xs, ys = ring.coords.xy
----> 7 xs.append(xs[1])
      8 ys.append(ys[1])
      9 return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(ring.coords)))/2.0

IndexError: array index out of range

:Overlay
   .WMTS.I  :WMTS   [Longitude,Latitude]
   .Image.I :Image   [x,y]   (value)

Summary

Starting from raw Landsat data, we have used a simple spectral clustering approach to observe the change in a lake water’s extent across time.

What’s next?

Adapt this notebook for your own use case or select another workflow example notebook.

Resources and References

Authored by Demetris Roumis circa Jan, 2023

Spectral Clustering

Overview

Prerequisites

Imports

Loading Data

Search the catalog

Load a dataset

Retain CRS Attribute

Reshaping Data

Standardize Data

ML pipeline

Client

Cluster Info

LocalCluster

Scheduler Info

Scheduler

Workers

Worker: 0

Un-flattening

Spectral Clustering for 1988

Load the data

Reshape and Standardize

Spectral Clustering

Unstack and Visualize

Spectral Clustering Over Time

Summary

What’s next?

Resources and References