Spectral Clustering

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 - Planetary Computer

Necessary

scikit-learn

Helpful

Spectral clustering

dask_ml

Helpful

Spectral clustering at scale

  • Time to learn: 20 minutes.

Imports

# Data
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

# Analysis
from dask_ml.cluster import SpectralClustering

# 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-6060ada8-f06c-11ee-8a9f-6045bdc85b55

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

Cluster Info

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:3157: 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 22 s, sys: 13.2 s, total: 35.2 s
Wall time: 29.2 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.

Un-flattening

Once the computation is done, we can use the coordinates of our input array to restack our output array back into an image. Again, one of the main benefits of using xarray for this stacking and unstacking is that it keeps track of the coordinate information for us.

Since the original array is n_samples by n_features (90000, 6) and the cluster label output is (90000,), we just need the coordinates from one of the original features in the shape of n_samples. We can just copy the coordinates from the first input feature and populate is with our output data:

template = flattened_t_xda[:, 0]
output_array = template.copy(data=labels)
output_array
<xarray.DataArray (z: 1071600)>
array([0, 2, 2, ..., 2, 2, 2], dtype=int32)
Coordinates:
    spatial_ref  int32 32611
    time         datetime64[ns] 2017-07-02T18:33:06.200763
    band         <U3 'red'
  * 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

With this new output array with coordinates copied from the input array, we can unstack back to the original x and y image dimensions by just using .unstack().

unstacked_2017 = output_array.unstack()
unstacked_2017
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/numpy/core/numeric.py:407: RuntimeWarning: invalid value encountered in cast
  multiarray.copyto(res, fill_value, casting='unsafe')
<xarray.DataArray (x: 950, y: 1128)>
array([[0, 2, 2, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [2, 0, 2, ..., 0, 0, 0],
       ...,
       [2, 2, 2, ..., 2, 2, 2],
       [2, 2, 2, ..., 2, 2, 2],
       [2, 2, 0, ..., 2, 2, 2]], dtype=int32)
Coordinates:
  * x            (x) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.638e+05
  * y            (y) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
    spatial_ref  int32 32611
    time         datetime64[ns] 2017-07-02T18:33:06.200763
    band         <U3 'red'
Attributes:
    crs:      epsg:32611

Finally, we can visualize the results! By hovering over the resulting imge, we can see that the lake water has been clustered with a certain label or ‘value’.

raw_plot_2017 = da_2017.sel(band="red").hvplot.image(
    x="x", y="y", geo=True, xlabel="lon", ylabel="lat", datashade=True, cmap="greys", title="Raw Image 2017",
)

result_plot_2017 = unstacked_2017.hvplot(
    x="x", y="y", cmap="Set3", geo=True, xlabel="lon", ylabel="lat", colorbar=False,  title="Spectral Clustering 2017",
)

raw_plot_2017 + result_plot_2017
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/operation/__init__.py:14: HoloviewsDeprecationWarning: 'ResamplingOperation' is deprecated and will be removed in version 1.18, use 'ResampleOperation2D' instead.
  from holoviews.operation.datashader import (

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:3157: 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.4 s, sys: 2.57 s, total: 21.9 s
Wall time: 26.8 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)

Unstack and Visualize

template = flattened_t_xda[:, 0]
output_array = template.copy(data=labels)
unstacked_1988 = output_array.unstack()
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/numpy/core/numeric.py:407: RuntimeWarning: invalid value encountered in cast
  multiarray.copyto(res, fill_value, casting='unsafe')
unstacked_1988
<xarray.DataArray (x: 950, y: 1128)>
array([[2, 0, 0, ..., 0, 0, 0],
       [2, 2, 2, ..., 0, 0, 0],
       [2, 2, 2, ..., 0, 0, 0],
       ...,
       [0, 0, 0, ..., 2, 2, 2],
       [0, 0, 0, ..., 2, 2, 2],
       [0, 0, 0, ..., 2, 2, 2]], dtype=int32)
Coordinates:
  * x            (x) float64 3.353e+05 3.353e+05 ... 3.637e+05 3.638e+05
  * y            (y) float64 4.301e+06 4.301e+06 ... 4.267e+06 4.267e+06
    spatial_ref  int32 32611
    time         datetime64[ns] 1988-07-02T18:03:59.010013
    band         <U3 'red'
Attributes:
    crs:      epsg:32611
raw_plot_1988 = da_1988.sel(band="red").hvplot.image(
    x="x", y="y", geo=True, xlabel="lon", ylabel="lat", datashade=True, cmap="greys", title="Raw 1988"
)

result_plot_1988 = unstacked_1988.hvplot(
    x="x", y="y", cmap="Set3", geo=True, xlabel="lon", ylabel="lat", colorbar=False, title="Spectral Clustering 1988",
)

raw_plot_1988 + result_plot_1988

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

By hovering over the lake in the center of each image, we can see that the water was labeled with a cluster label value in both images. Let’s programmatically grab the cluster label at the center for each image.

def get_center_value(arr):
    center_y = arr.shape[0] // 2
    center_x = arr.shape[1] // 2
    center_value = arr[center_y, center_x]
    return int(center_value.values)

water_cluster_1988_label = get_center_value(result_plot_1988.data.value)
water_cluster_2017_label = get_center_value(result_plot_2017.data.value)

Now, we want to align the cluster label for water. Let’s try to set anything that is water to 1 and otherwise 0.

water_1988 = xr.where(unstacked_1988 == water_cluster_1988_label, 1, 0)
water_2017 = xr.where(unstacked_2017 == water_cluster_2017_label, 1, 0)
water_1988.hvplot(geo=True)
---------------------------------------------------------------------------
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:1287, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1280 def _repr_mimebundle_(self, include=None, exclude=None):
   1281     """
   1282     Resolves the class hierarchy for the class rendering the
   1283     object using any display hooks registered on Store.display
   1284     hooks.  The output of all registered display_hooks is then
   1285     combined and returned.
   1286     """
-> 1287     return Store.render(self)

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

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

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

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

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

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

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

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:509, in Renderable._render_model(self, doc, comm)
    507 if comm is None:
    508     comm = state._comm_manager.get_server_comm()
--> 509 model = self.get_root(doc, comm)
    511 if config.embed:
    512     embed_state(self, model, doc,
    513                 json=config.embed_json,
    514                 json_prefix=config.embed_json_prefix,
    515                 save_path=config.embed_save_path,
    516                 load_path=config.embed_load_path,
    517                 progress=False)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:560, in Renderable.get_root(self, doc, comm, preprocess)
    543 """
    544 Returns the root model and applies pre-processing hooks
    545 
   (...)
    557 Returns the bokeh model corresponding to this panel object
    558 """
    559 doc = init_doc(doc)
--> 560 root = self._get_model(doc, comm=comm)
    561 if preprocess:
    562     self._preprocess(root)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:146, in Panel._get_model(self, doc, root, parent, comm)
    144 if root is None:
    145     root = model
--> 146 objects = self._get_objects(model, [], doc, root, comm)
    147 props = dict(self._init_params(), objects=objects)
    148 model.update(**self._process_param_change(props))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:131, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    129 else:
    130     try:
--> 131         child = pane._get_model(doc, root, model, comm)
    132     except RerenderError:
    133         return self._get_objects(model, current_objects[:i], doc, root, comm)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:367, in HoloViews._get_model(self, doc, root, parent, comm)
    365     plot = self.object
    366 else:
--> 367     plot = self._render(doc, comm, root)
    369 plot.pane = self
    370 backend = plot.renderer.backend

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/renderer.py:70, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     63 @bothmethod
     64 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     65     """
     66     Given a HoloViews Viewable return a corresponding plot instance.
     67     Allows supplying a document attach the plot to, useful when
     68     combining the bokeh model with another plot.
     69     """
---> 70     plot = super().get_plot(obj, doc, renderer, **kwargs)
     71     if plot.document is None:
     72         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:241, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    238     defaults = [kd.default for kd in plot.dimensions]
    239     init_key = tuple(v if d is None else d for v, d in
    240                      zip(plot.keys[0], defaults))
--> 241     plot.update(init_key)
    242 else:
    243     plot = obj

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

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

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

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:426, in ElementPlot._axis_props(self, plots, subplots, element, ranges, pos, dim, range_tags_extras, extra_range_name)
    424     l, b, r, t = self.get_extents(range_el, ranges, dim)
    425 else:
--> 426     l, b, r, t = self.get_extents(range_el, ranges)
    427 if self.invert_axes:
    428     l, b, r, t = b, l, t, r

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/plot.py:73, in ProjectionPlot.get_extents(self, element, ranges, range_type)
     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:103, in project_extents(extents, src_proj, dest_proj, tol)
    101     geom_in_src_proj = geom_clipped_to_dest_proj
    102 try:
--> 103     geom_in_crs = dest_proj.project_geometry(geom_in_src_proj, src_proj)
    104 except ValueError:
    105     src_name =type(src_proj).__name__

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:953, in Projection._project_polygon(self, polygon, src_crs)
    951     is_ccw = True
    952 else:
--> 953     is_ccw = polygon.exterior.is_ccw
    954 # Project the polygon exterior/interior rings.
    955 # Each source ring will result in either a ring, or one or more
    956 # lines.
    957 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   [y,x]   (value)
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:1287, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1280 def _repr_mimebundle_(self, include=None, exclude=None):
   1281     """
   1282     Resolves the class hierarchy for the class rendering the
   1283     object using any display hooks registered on Store.display
   1284     hooks.  The output of all registered display_hooks is then
   1285     combined and returned.
   1286     """
-> 1287     return Store.render(self)

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

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

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

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

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

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

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

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:509, in Renderable._render_model(self, doc, comm)
    507 if comm is None:
    508     comm = state._comm_manager.get_server_comm()
--> 509 model = self.get_root(doc, comm)
    511 if config.embed:
    512     embed_state(self, model, doc,
    513                 json=config.embed_json,
    514                 json_prefix=config.embed_json_prefix,
    515                 save_path=config.embed_save_path,
    516                 load_path=config.embed_load_path,
    517                 progress=False)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:560, in Renderable.get_root(self, doc, comm, preprocess)
    543 """
    544 Returns the root model and applies pre-processing hooks
    545 
   (...)
    557 Returns the bokeh model corresponding to this panel object
    558 """
    559 doc = init_doc(doc)
--> 560 root = self._get_model(doc, comm=comm)
    561 if preprocess:
    562     self._preprocess(root)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:146, in Panel._get_model(self, doc, root, parent, comm)
    144 if root is None:
    145     root = model
--> 146 objects = self._get_objects(model, [], doc, root, comm)
    147 props = dict(self._init_params(), objects=objects)
    148 model.update(**self._process_param_change(props))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:131, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    129 else:
    130     try:
--> 131         child = pane._get_model(doc, root, model, comm)
    132     except RerenderError:
    133         return self._get_objects(model, current_objects[:i], doc, root, comm)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:367, in HoloViews._get_model(self, doc, root, parent, comm)
    365     plot = self.object
    366 else:
--> 367     plot = self._render(doc, comm, root)
    369 plot.pane = self
    370 backend = plot.renderer.backend

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/renderer.py:70, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     63 @bothmethod
     64 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     65     """
     66     Given a HoloViews Viewable return a corresponding plot instance.
     67     Allows supplying a document attach the plot to, useful when
     68     combining the bokeh model with another plot.
     69     """
---> 70     plot = super().get_plot(obj, doc, renderer, **kwargs)
     71     if plot.document is None:
     72         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:241, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    238     defaults = [kd.default for kd in plot.dimensions]
    239     init_key = tuple(v if d is None else d for v, d in
    240                      zip(plot.keys[0], defaults))
--> 241     plot.update(init_key)
    242 else:
    243     plot = obj

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

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

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

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

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

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:426, in ElementPlot._axis_props(self, plots, subplots, element, ranges, pos, dim, range_tags_extras, extra_range_name)
    424     l, b, r, t = self.get_extents(range_el, ranges, dim)
    425 else:
--> 426     l, b, r, t = self.get_extents(range_el, ranges)
    427 if self.invert_axes:
    428     l, b, r, t = b, l, t, r

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/plot.py:73, in ProjectionPlot.get_extents(self, element, ranges, range_type)
     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:103, in project_extents(extents, src_proj, dest_proj, tol)
    101     geom_in_src_proj = geom_clipped_to_dest_proj
    102 try:
--> 103     geom_in_crs = dest_proj.project_geometry(geom_in_src_proj, src_proj)
    104 except ValueError:
    105     src_name =type(src_proj).__name__

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:953, in Projection._project_polygon(self, polygon, src_crs)
    951     is_ccw = True
    952 else:
--> 953     is_ccw = polygon.exterior.is_ccw
    954 # Project the polygon exterior/interior rings.
    955 # Each source ring will result in either a ring, or one or more
    956 # lines.
    957 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:1287, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1280 def _repr_mimebundle_(self, include=None, exclude=None):
   1281     """
   1282     Resolves the class hierarchy for the class rendering the
   1283     object using any display hooks registered on Store.display
   1284     hooks.  The output of all registered display_hooks is then
   1285     combined and returned.
   1286     """
-> 1287     return Store.render(self)

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

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

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

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

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

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

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

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:509, in Renderable._render_model(self, doc, comm)
    507 if comm is None:
    508     comm = state._comm_manager.get_server_comm()
--> 509 model = self.get_root(doc, comm)
    511 if config.embed:
    512     embed_state(self, model, doc,
    513                 json=config.embed_json,
    514                 json_prefix=config.embed_json_prefix,
    515                 save_path=config.embed_save_path,
    516                 load_path=config.embed_load_path,
    517                 progress=False)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:560, in Renderable.get_root(self, doc, comm, preprocess)
    543 """
    544 Returns the root model and applies pre-processing hooks
    545 
   (...)
    557 Returns the bokeh model corresponding to this panel object
    558 """
    559 doc = init_doc(doc)
--> 560 root = self._get_model(doc, comm=comm)
    561 if preprocess:
    562     self._preprocess(root)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:146, in Panel._get_model(self, doc, root, parent, comm)
    144 if root is None:
    145     root = model
--> 146 objects = self._get_objects(model, [], doc, root, comm)
    147 props = dict(self._init_params(), objects=objects)
    148 model.update(**self._process_param_change(props))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:131, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    129 else:
    130     try:
--> 131         child = pane._get_model(doc, root, model, comm)
    132     except RerenderError:
    133         return self._get_objects(model, current_objects[:i], doc, root, comm)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:367, in HoloViews._get_model(self, doc, root, parent, comm)
    365     plot = self.object
    366 else:
--> 367     plot = self._render(doc, comm, root)
    369 plot.pane = self
    370 backend = plot.renderer.backend

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/renderer.py:70, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     63 @bothmethod
     64 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     65     """
     66     Given a HoloViews Viewable return a corresponding plot instance.
     67     Allows supplying a document attach the plot to, useful when
     68     combining the bokeh model with another plot.
     69     """
---> 70     plot = super().get_plot(obj, doc, renderer, **kwargs)
     71     if plot.document is None:
     72         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:241, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    238     defaults = [kd.default for kd in plot.dimensions]
    239     init_key = tuple(v if d is None else d for v, d in
    240                      zip(plot.keys[0], defaults))
--> 241     plot.update(init_key)
    242 else:
    243     plot = obj

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

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

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

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:426, in ElementPlot._axis_props(self, plots, subplots, element, ranges, pos, dim, range_tags_extras, extra_range_name)
    424     l, b, r, t = self.get_extents(range_el, ranges, dim)
    425 else:
--> 426     l, b, r, t = self.get_extents(range_el, ranges)
    427 if self.invert_axes:
    428     l, b, r, t = b, l, t, r

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geoviews/plotting/plot.py:73, in ProjectionPlot.get_extents(self, element, ranges, range_type)
     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:103, in project_extents(extents, src_proj, dest_proj, tol)
    101     geom_in_src_proj = geom_clipped_to_dest_proj
    102 try:
--> 103     geom_in_crs = dest_proj.project_geometry(geom_in_src_proj, src_proj)
    104 except ValueError:
    105     src_name =type(src_proj).__name__

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:953, in Projection._project_polygon(self, polygon, src_crs)
    951     is_ccw = True
    952 else:
--> 953     is_ccw = polygon.exterior.is_ccw
    954 # Project the polygon exterior/interior rings.
    955 # Each source ring will result in either a ring, or one or more
    956 # lines.
    957 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:1287, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1280 def _repr_mimebundle_(self, include=None, exclude=None):
   1281     """
   1282     Resolves the class hierarchy for the class rendering the
   1283     object using any display hooks registered on Store.display
   1284     hooks.  The output of all registered display_hooks is then
   1285     combined and returned.
   1286     """
-> 1287     return Store.render(self)

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

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

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

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

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

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

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

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:509, in Renderable._render_model(self, doc, comm)
    507 if comm is None:
    508     comm = state._comm_manager.get_server_comm()
--> 509 model = self.get_root(doc, comm)
    511 if config.embed:
    512     embed_state(self, model, doc,
    513                 json=config.embed_json,
    514                 json_prefix=config.embed_json_prefix,
    515                 save_path=config.embed_save_path,
    516                 load_path=config.embed_load_path,
    517                 progress=False)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/viewable.py:560, in Renderable.get_root(self, doc, comm, preprocess)
    543 """
    544 Returns the root model and applies pre-processing hooks
    545 
   (...)
    557 Returns the bokeh model corresponding to this panel object
    558 """
    559 doc = init_doc(doc)
--> 560 root = self._get_model(doc, comm=comm)
    561 if preprocess:
    562     self._preprocess(root)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:146, in Panel._get_model(self, doc, root, parent, comm)
    144 if root is None:
    145     root = model
--> 146 objects = self._get_objects(model, [], doc, root, comm)
    147 props = dict(self._init_params(), objects=objects)
    148 model.update(**self._process_param_change(props))

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/layout/base.py:131, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    129 else:
    130     try:
--> 131         child = pane._get_model(doc, root, model, comm)
    132     except RerenderError:
    133         return self._get_objects(model, current_objects[:i], doc, root, comm)

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/panel/pane/holoviews.py:367, in HoloViews._get_model(self, doc, root, parent, comm)
    365     plot = self.object
    366 else:
--> 367     plot = self._render(doc, comm, root)
    369 plot.pane = self
    370 backend = plot.renderer.backend

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/renderer.py:70, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     63 @bothmethod
     64 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     65     """
     66     Given a HoloViews Viewable return a corresponding plot instance.
     67     Allows supplying a document attach the plot to, useful when
     68     combining the bokeh model with another plot.
     69     """
---> 70     plot = super().get_plot(obj, doc, renderer, **kwargs)
     71     if plot.document is None:
     72         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:241, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    238     defaults = [kd.default for kd in plot.dimensions]
    239     init_key = tuple(v if d is None else d for v, d in
    240                      zip(plot.keys[0], defaults))
--> 241     plot.update(init_key)
    242 else:
    243     plot = obj

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

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

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

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/holoviews/plotting/bokeh/element.py:424, in ElementPlot._axis_props(self, plots, subplots, element, ranges, pos, dim, range_tags_extras, extra_range_name)
    422 else:
    423     if isinstance(self, OverlayPlot):
--> 424         l, b, r, t = self.get_extents(range_el, ranges, dim)
    425     else:
    426         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)
     65     (x0, x1), (y0, y1) = proj.x_limits, proj.y_limits
     66     return (x0, y0, x1, y1)
---> 67 extents = super(ProjectionPlot, self).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:1970, in GenericOverlayPlot.get_extents(self, overlay, ranges, range_type, dimension, **kwargs)
   1969 def get_extents(self, overlay, ranges, range_type='combined', dimension=None, **kwargs):
-> 1970     subplot_extents = self._get_subplot_extents(overlay, ranges, range_type, dimension)
   1971     zrange = isinstance(self.projection, str) and self.projection == '3d'
   1972     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:1965, in GenericOverlayPlot._get_subplot_extents(self, overlay, ranges, range_type, dimension)
   1963         sp_ranges = util.match_spec(layer, ranges) if ranges else {}
   1964     for rt in extents:
-> 1965         extent = subplot.get_extents(layer, sp_ranges, range_type=rt)
   1966         extents[rt].append(extent)
   1967 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)
     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:103, in project_extents(extents, src_proj, dest_proj, tol)
    101     geom_in_src_proj = geom_clipped_to_dest_proj
    102 try:
--> 103     geom_in_crs = dest_proj.project_geometry(geom_in_src_proj, src_proj)
    104 except ValueError:
    105     src_name =type(src_proj).__name__

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

File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/crs.py:953, in Projection._project_polygon(self, polygon, src_crs)
    951     is_ccw = True
    952 else:
--> 953     is_ccw = polygon.exterior.is_ccw
    954 # Project the polygon exterior/interior rings.
    955 # Each source ring will result in either a ring, or one or more
    956 # lines.
    957 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