Data Ingestion - Geospatial-Specific Tooling

PySTAC

Overview

In this notebook, you will ingest Landsat data for use in machine learning. Machine learning tasks often involve a lot of data, and in Python, data is typically stored in memory as simple NumPy arrays. However, higher-level containers built on top of NumPy arrays provide more functionality for multidimensional gridded data (xarray) or out-of-core and distributed data (Dask). Our goal for data ingestion will be to load specific Landsat data of interest into one of these higher-level containers.

Microsoft Plantery Computer is one of several providers of Landsat Data. We are using it together with pystac-client and odc-stac because together they provide a nice Python API for searching and loading with specific criteria such as spatial area, datetime, Landsat mission, and cloud coverage.

Earth science datasets are often stored on remote servers that may be too large to download locally. Therefore, in this cookbook, we will focus primarily on ingestion approaches that load small portions of data from a remote source, as needed. However, the approach for your own work will depend not only on data size and location but also the intended analysis, so in a follow up notebook, you will see an alternative approache for generalized data access and management.

Prerequisites

Concepts

Importance

Notes

Intro to Landsat

Necessary

Background

About the Microsoft Planetary Computer

Helpful

Background

pystac-client Usage

Helpful

Consult as needed

odc.stac.load Reference

Helpful

Consult as needed

xarray

Necessary

Intro to Dask Array

Helpful

Panel Getting Started Guide

Helpful

  • Time to learn: 10 minutes

Imports

import odc.stac
import pandas as pd
import planetary_computer
import pystac_client
import xarray as xr
from pystac.extensions.eo import EOExtension as eo

# Viz
import hvplot.xarray
import panel as pn

pn.extension()

Open and read the root of the STAC catalog

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

'Microsoft Planetary Computer STAC API'

Microsoft Planetary Computer has a public STAC metadata but the actual data assets are in private Azure Blob Storage containers and require authentication. pystac-client provides a modifier keyword that we can use to manually sign the item. Otherwise, we’d get an error when trying to access the asset.

Search for Landsat Data

Let’s say that an analysis we want to run requires landsat data over a specific region and from a specific time period. We can use our catalog to search for assets that fit our search criteria.

First, let’s find the name of the landsat dataset. This page is a nice resource for browsing the available collections, but we can also just search the catalog for ‘landsat’:

all_collections = [i.id for i in catalog.get_collections()]
landsat_collections = [
    collection for collection in all_collections if "landsat" in collection
]
landsat_collections

['landsat-c2-l2', 'landsat-c2-l1']

We’ll use the landsat-c2-l2 dataset, which stands for Collection 2 Level-2. It contains data from several landsat missions and has better data quality than Level 1 (landsat-c2-l1). Microsoft Planetary Computer has descriptions of Level 1 and Level 2, but a direct and succinct comparison can be found in this community post, and the information can be verified with USGS.

Now, let’s set our search parameters. You may already know the bounding box (region/area of interest) coordinates, but if you don’t, there are many useful tools like bboxfinder.com that can help.

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"

We can also specify other parameters in the query, such as a specific landsat mission and the max percent of cloud cover:

platform = "landsat-8"
cloudy_less_than = 1  # percent

Now we run the search and list the results:

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

Returned 3 Items:

{(0, 'LC08_L2SP_042033_20170718_02_T1'): 0,
 (1, 'LC08_L2SP_042033_20170702_02_T1'): 1,
 (2, 'LC08_L2SP_042033_20170616_02_T1'): 2}

It looks like there were three image stacks taken by Landsat 8 over this spatial region during the summer months of 2017 that has less than 1 percent cloud cover.

Preview Results and Select a Dataset

Before loading one of the available image stacks, it would be useful to get a visual check of the results. Many datasets have a rendered preview or thumbnail image that can be accessed without having to load the full resolution data.

We can create a simple interactive application using the Panel library to access and display rendered PNG previews of the our search results. Note that these pre-rendered images are of large tiles that span beyond our bounding box of interest. In the next steps, we will only be loading in a small area around the lake.

item_sel = pn.widgets.Select(value=1, options=item_id, name="item")

def get_preview(i):
    return pn.panel(items[i].assets["rendered_preview"].href, height=300)


pn.Row(item_sel, pn.bind(get_preview, item_sel))

selected_item = items[1]
selected_item

Access the Data

Now that we have selected a dataset from our catalog, we can procede to access the data. We want to be very selective about the data that we read and when we read it because the amount of downloaded data can quickly get out of hand. Therefore, let’s select only a subset of images.

First, we’ll preview the different image assets (or Bands) available in the Landsat item.

assets = []
for _, asset in selected_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

Then we will select a few bands (images) of interest:

bands_of_interest = ["red", "green", "blue"]

Finally, we lazily load the selected data. We will use the package called odc which allows us to load only a specific region of interest (bounding box or ‘bbox’) and specific bands (images) of interest. We will also use the chunks argument to load the data as dask arrays; this will load the metadata now and delay the loading until we actually use the data, or until we force the data to be loaded by using .compute().

ds = odc.stac.stac_load(
    [selected_item],
    bands=bands_of_interest,
    bbox=bbox,
    chunks={},  # <-- use Dask
).isel(time=0)
ds

<xarray.Dataset>
Dimensions:      (y: 1128, x: 950)
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
Data variables:
    red          (y, x) uint16 dask.array<chunksize=(1128, 950), meta=np.ndarray>
    green        (y, x) uint16 dask.array<chunksize=(1128, 950), meta=np.ndarray>
    blue         (y, x) uint16 dask.array<chunksize=(1128, 950), meta=np.ndarray>

Let’s combine the bands of the dataset into a single DataArray that has the band names as coordinates of a new ‘band’ dimension, and also call .compute() to finally load the data.

da = ds.to_array(dim="band").compute()
da

<xarray.DataArray (band: 3, y: 1128, x: 950)>
array([[[14691, 14914, 14988, ..., 16283, 16292, 16316],
        [14655, 14859, 14969, ..., 16272, 16185, 16079],
        [14531, 14699, 14972, ..., 15318, 15526, 14734],
        ...,
        [13804, 13561, 13601, ..., 18311, 18202, 17625],
        [13857, 13828, 13858, ..., 19400, 18942, 18551],
        [13840, 13786, 13867, ..., 17873, 17917, 18453]],

       [[13233, 13402, 13565, ..., 14553, 14658, 14657],
        [13291, 13428, 13585, ..., 14590, 14478, 14550],
        [13122, 13287, 13601, ..., 13987, 14220, 13571],
        ...,
        [12720, 12552, 12468, ..., 16580, 16411, 15899],
        [12704, 12644, 12658, ..., 17351, 16853, 16505],
        [12647, 12620, 12698, ..., 15990, 16211, 16686]],

       [[11572, 11629, 11723, ..., 12857, 12918, 12946],
        [11588, 11655, 11721, ..., 12848, 12792, 12715],
        [11510, 11608, 11781, ..., 12371, 12453, 12053],
        ...,
        [11195, 11104, 11045, ..., 14182, 14031, 13716],
        [11125, 11061, 11106, ..., 14652, 14284, 14062],
        [11059, 11050, 11134, ..., 13756, 13865, 14209]]], dtype=uint16)
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' 'green' 'blue'

Visualize the data

Often, data ingestion involves quickly visualizing your raw data to get a sense that things are proceeding accordingly. As we have created an array with red, blue, and green bands, we can quickly display a natural color image of the lake using the .plot.imshow() function of xarray. We’ll use the robust=True argument because the data values are outside the range of typical RGB images.

da.plot.imshow(robust=True, size=3)

<matplotlib.image.AxesImage at 0x7f5b18271000>
../_images/44752ad5ede4bfd4ff006172839266a1d5ca899c87f3d555f5c61208891dfbd7.png

Now, let’s use hvplot to provide an interactive visualization of the inividual bands in our array.

ds

<xarray.Dataset>
Dimensions:      (y: 1128, x: 950)
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
Data variables:
    red          (y, x) uint16 dask.array<chunksize=(1128, 950), meta=np.ndarray>
    green        (y, x) uint16 dask.array<chunksize=(1128, 950), meta=np.ndarray>
    blue         (y, x) uint16 dask.array<chunksize=(1128, 950), meta=np.ndarray>
da.hvplot.image(x="x", y="y", cmap="viridis", aspect=1)

Let’s plot the bands as seperate columns by specifying a dimension to expand with col='band'. We can also set rasterize=True to use Datashader (another HoloViz tool) to render large data into a 2D histogram, where every array cell counts the data points falling into that pixel, as set by the resolution of your screen. This is especially important for large and high resolution images that would otherwise cause issues when attempting to render in a browser.

da.hvplot.image(
    x="x", y="y", col="band", cmap="viridis", xaxis=False, yaxis=False, colorbar=False, rasterize=True
)

/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(
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/geopandas/_compat.py:154: UserWarning: The Shapely GEOS version (3.11.1-CAPI-1.17.1) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow.
  set_use_pygeos()

Select the zoom tool and zoom in on of the plots to see that all the images are all automatically linked!

Retain Attributes

When working with many image arrays, it’s critical to retain the data properties as xarray attributes:

da.attrs = selected_item.properties
da

<xarray.DataArray (band: 3, y: 1128, x: 950)>
array([[[14691, 14914, 14988, ..., 16283, 16292, 16316],
        [14655, 14859, 14969, ..., 16272, 16185, 16079],
        [14531, 14699, 14972, ..., 15318, 15526, 14734],
        ...,
        [13804, 13561, 13601, ..., 18311, 18202, 17625],
        [13857, 13828, 13858, ..., 19400, 18942, 18551],
        [13840, 13786, 13867, ..., 17873, 17917, 18453]],

       [[13233, 13402, 13565, ..., 14553, 14658, 14657],
        [13291, 13428, 13585, ..., 14590, 14478, 14550],
        [13122, 13287, 13601, ..., 13987, 14220, 13571],
        ...,
        [12720, 12552, 12468, ..., 16580, 16411, 15899],
        [12704, 12644, 12658, ..., 17351, 16853, 16505],
        [12647, 12620, 12698, ..., 15990, 16211, 16686]],

       [[11572, 11629, 11723, ..., 12857, 12918, 12946],
        [11588, 11655, 11721, ..., 12848, 12792, 12715],
        [11510, 11608, 11781, ..., 12371, 12453, 12053],
        ...,
        [11195, 11104, 11045, ..., 14182, 14031, 13716],
        [11125, 11061, 11106, ..., 14652, 14284, 14062],
        [11059, 11050, 11134, ..., 13756, 13865, 14209]]], dtype=uint16)
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' 'green' 'blue'
Attributes: (12/22)
    gsd:                          30
    created:                      2022-05-06T17:46:34.110946Z
    sci:doi:                      10.5066/P9OGBGM6
    datetime:                     2017-07-02T18:33:06.200763Z
    platform:                     landsat-8
    proj:epsg:                    32611
    ...                           ...
    view:sun_azimuth:             125.03739105
    landsat:correction:           L2SP
    view:sun_elevation:           65.85380157
    landsat:cloud_cover_land:     0.53
    landsat:collection_number:    02
    landsat:collection_category:  T1

Notice that you can now expand the Attributes: dropdown to see the properties of this data.

Set the crs attribute

As the data is in ‘meter’ units from a reference point, we can plot in commonly used longitude, latitude coordinates with .hvplot(geo=True) if our array has a valid coordinate reference system (CRS) attribute. This value is provided from Microsoft Planetary Computer as the proj:epsg property, so we just need to copy it to a new attribute crs so that hvPlot can automatically find it, without us having to further specify anything in our plotting code

Note, this CRS is referenced by an EPSG code that can be accessed from the metadata of our selected catalog search result. We can see more about this dataset’s specific code at EPSG.io/32611. You can also read more about EPSG codes in general in this Coordinate Reference Systems: EPSG codes online book chapter.

da.attrs["crs"] = f"epsg:{selected_item.properties['proj:epsg']}"
da.attrs["crs"]

'epsg:32611'

Now we can use .hvplot(geo=True) to plot in longitude and latitude coordinates. Informing hvPlot that this is geographic data also allows us to overlay data on aligned geographic tiles using the tiles parameter.

da.hvplot.image(
    x="x", y="y", cmap="viridis", geo=True, alpha=.9, tiles="ESRI", xlabel="Longitude", ylabel="Latitude", colorbar=False, aspect=1,
)

---------------------------------------------------------------------------
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:261, in display(obj, raw_output, **kwargs)
    259 elif isinstance(obj, (HoloMap, DynamicMap)):
    260     with option_state(obj):
--> 261         output = map_display(obj)
    262 elif isinstance(obj, Plot):
    263     output = render(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:209, in map_display(vmap, max_frames)
    206     max_frame_warning(max_frames)
    207     return None
--> 209 return render(vmap)

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

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

Summary

The data access approach should adapt to features of the data and your intended analysis. As Landsat data is large and multidimensional, a good approach is to use Microsoft Plantery Computer, pystac-client, and odc-stac together for searching the metadata catalog and lazily loading specific data chunks. Once you have accessed data, visualize it with hvPlot to ensure that it matches your expectations.

What’s next?

Before we proceed to workflow examples, we can explore an alternate way of accessing data using generalized tooling.

Resources and References

  • Authored by Demetris Roumis circa Jan, 2023

  • Guidance for parts of this notebook was provided by Microsoft in ‘Reading Data from the STAC API’

  • The image used in the banner is from an announcement about PySTAC from Azavea

On this page