Subset 3D data using rectangle tool and plot it with time control

Hi All,

I have a 3D array (lat, lon, time) and want to achieve the function like this:

  1. Show the World map without any data shown
  2. Enlarge and use the rectangle tool to select one region
  3. Load the data by the selected region and plot the first-time data.
  4. Using the time control bar to see the data changes.

Here’s the example data

import xarray as xr
air_ds = xr.tutorial.open_dataset('air_temperature').load()
print(air_ds)
<xarray.Dataset>
Dimensions:  (lat: 25, lon: 53, time: 2920)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Data variables:
    air      (time, lat, lon) float32 241.2 242.5 243.5 ... 296.5 296.2 295.7

I would appreciate it a lot if you can give me any advice/examples
Thanks!

Here’s how I coded something like this
Step 1 import + make the base map (1.)

import panel as pn
import xarray as xr
import holoviews as hv

hv.extension("bokeh")

tiles = hv.element.CartoDark()

Step 2: instantiate a stream (2.)

import panel as pn
import xarray as xr
import holoviews as hv

hv.extension("bokeh")

tiles = hv.element.CartoDark()
stream = hv.streams.SelectionXY()

Step 3: do something with the stream and add a basic print

import panel as pn
import xarray as xr
import holoviews as hv

hv.extension("bokeh")

def plot(x):
    print(x)

tiles = hv.element.CartoDark()
stream = hv.streams.SelectionXY()
dmap = hv.DynamicMap(plot, streams=[stream])

tiles * dmap

Step 4: keep crashing and revising
KeyError: "Callable 'plot' missing keywords to accept stream parameters: bounds,x_selection,y_selection"
TypeError: DynamicMap does not accept NoneType type, data elements have to be a ('ViewableElement', 'NdMapping', 'Layout').

import panel as pn
import xarray as xr
import holoviews as hv

hv.extension("bokeh")

def plot(bounds, x_selection, y_selection):
    print(bounds, x_selection, y_selection)
    return hv.Image([])

tiles = hv.element.CartoDark()
stream = hv.streams.SelectionXY()
dmap = hv.DynamicMap(plot, streams=[stream])

tiles * dmap

Step 5: test out stream input (zoom out + use the selection tool).

Step 6: realize that I don’t want to deal with projecting Mercator northing/easting to degrees for an example :smiley:; swap to basic coastlines projected in degrees and test inputs again

import panel as pn
import xarray as xr
import holoviews as hv
import geoviews as gv
import cartopy.crs as ccrs

gv.extension("bokeh")

def plot(bounds, x_selection, y_selection):
    print(bounds, x_selection, y_selection)
    return hv.Image([])

coastline = gv.feature.coastline().opts(
    projection=ccrs.PlateCarree(), global_extent=True)
stream = hv.streams.SelectionXY()
dmap = gv.DynamicMap(plot, streams=[stream])

coastline * dmap

Step 7: now we’re ready to fill in the function (3.)

import panel as pn
import xarray as xr
import holoviews as hv
import geoviews as gv
import cartopy.crs as ccrs

gv.extension("bokeh")
da = xr.tutorial.open_dataset("air_temperature")["air"].isel(time=0)
da["lon"] = (da["lon"] + 180) % 360 - 180
da = da.sortby("lat")

def plot(bounds, x_selection, y_selection):
    if x_selection is not None:
        da_sel = da.sel(lon=slice(*x_selection), lat=slice(*y_selection))
        return hv.Image(da_sel, ["lon", "lat"])
    else:
        return hv.Image([])

coastline = gv.feature.coastline().opts(
    projection=ccrs.PlateCarree(), global_extent=True)
stream = hv.streams.SelectionXY()
dmap = gv.DynamicMap(plot, streams=[stream])

coastline * dmap

Tada!

2 Likes

Oh and finally step 4:

import panel as pn
import xarray as xr
import holoviews as hv
import geoviews as gv
import cartopy.crs as ccrs

gv.extension("bokeh")
da = xr.tutorial.open_dataset("air_temperature")["air"]
da["lon"] = (da["lon"] + 180) % 360 - 180
da = da.sortby("lat")

slider = pn.widgets.DiscreteSlider(options=da["time"].dt.strftime("%Y-%m-%d %H%MZ").values.tolist())

def plot(bounds, x_selection, y_selection, value):
    if x_selection is not None:
        da_sel = da.sel(lon=slice(*x_selection), lat=slice(*y_selection), time=value)
        return hv.Dataset(da_sel, kdims=["lon", "lat"]).to.image()
    else:
        return hv.Image([])

coastline = gv.feature.coastline().opts(
    projection=ccrs.PlateCarree(), global_extent=True, active_tools=["box_select"])
stream = hv.streams.SelectionXY()
dmap = gv.DynamicMap(plot, streams=[stream, slider.param.value])
overlay = coastline * dmap
pn.Row(overlay, slider)

ezgif-2-6590c5013ed3

1 Like

@ahuang11 Wow, thanks! That’s quite nice! I will play with my data using different tiles :wink:
BTW, is it possible to support using the mouse to zoom in and out?

1 Like

Yes, manually clicking the toolbar, or:
active_tools=["box_select", "wheel_zoom"]

@ahuang11 I tested the final code. No matter how I draw the rectangle, no data is shown …

Version info:

holoviews: 1.14.4
geoviews: 1.9.1
panel: 0.11.3

Update: If I use jupyter-notebook instead of jupyter-lab, then it works well.
Solution for Jupyter-lab: jupyter labextension install @pyviz/jupyterlab_pyviz

1 Like

@ahuang11 The CartoDark sounds better for some cases, is it possible to make it run successfully?

1 Like

Yes, just properly convert the northings/eastings to lat/lon before select. You can use cartopy or holoviews:
https://holoviews.org/reference_manual/holoviews.util.html#holoviews.util.transform.easting_northing_to_lon_lat

I tried this code:

import panel as pn
import xarray as xr
import holoviews as hv
import geoviews as gv
import cartopy.crs as ccrs

gv.extension("bokeh")
da = xr.tutorial.open_dataset("air_temperature")["air"]

da["lon"] = (da["lon"] + 180) % 360 - 180
da = da.sortby("lat")

slider = pn.widgets.DiscreteSlider(options=da["time"].dt.strftime("%Y-%m-%d %H%MZ").values.tolist())

def plot(bounds, x_selection, y_selection, value):
    if x_selection is not None:
        da_sel = da.sel(lon=slice(*x_selection), lat=slice(*y_selection), time=value)
        return hv.Dataset(da_sel, kdims=["lon", "lat"]).to.image()
    else:
        return hv.Image([])

coastline = gv.tile_sources.CartoDark(projection=ccrs.PlateCarree(), global_extent=True, active_tools=["box_select"])

stream = hv.streams.SelectionXY()
dmap = gv.DynamicMap(plot, streams=[stream, slider.param.value])
overlay = coastline * dmap
pn.Row(overlay, slider)

The data is shown, but the lon/lat is wrong and map can’t be are shown even though I haven’t select the region yet:

You can’t use projection=ccrs.PlateCarree() for a tile source I think? Not to mention, you need to put those in .opts() (I assume you’re getting a bunch of param warnings)?

If I change it to

coastline = gv.tile_sources.CartoDark(active_tools=["box_select"])

I got this image which has no map shown:

If changed to this one:

coastline = gv.tile_sources.CartoDark(active_tools=["box_select"]).opts(projection=ccrs.PlateCarree())

Similar thing:

Then, adding the selection leads to the same error:

Traceback (most recent call last):
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/backends/netCDF4_.py", line 104, in _getitem
    array = getitem(original_array, key)
  File "src/netCDF4/_netCDF4.pyx", line 4397, in netCDF4._netCDF4.Variable.__getitem__
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/netCDF4/utils.py", line 467, in _out_array_shape
    c = count[..., i].ravel()[0] # All elements should be identical.
IndexError: index 0 is out of bounds for axis 0 with size 0

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/pyviz_comms/__init__.py", line 316, in _handle_msg
    self._on_msg(msg)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/bokeh/callbacks.py", line 169, in on_msg
    raise e
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/bokeh/callbacks.py", line 161, in on_msg
    Stream.trigger(streams)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/streams.py", line 186, in trigger
    subscriber(**dict(union))
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/plot.py", line 244, in refresh
    raise e
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/plot.py", line 240, in refresh
    self._trigger_refresh(stream_key)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/plot.py", line 257, in _trigger_refresh
    self.update(key)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/plot.py", line 983, in update
    item = self.__getitem__(key)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/plot.py", line 446, in __getitem__
    self.update_frame(frame)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/bokeh/element.py", line 2449, in update_frame
    subplot.update_frame(key, ranges, element=el)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/bokeh/element.py", line 1538, in update_frame
    self._update_glyphs(element, ranges, self.style[self.cyclic_index])
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/bokeh/element.py", line 1445, in _update_glyphs
    data, mapping, style = self.get_data(element, ranges, style)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/plotting/bokeh/raster.py", line 112, in get_data
    img = element.dimension_values(i, flat=False)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/core/data/__init__.py", line 205, in pipelined_fn
    result = method_fn(*args, **kwargs)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/core/data/__init__.py", line 1105, in dimension_values
    values = self.interface.values(self, dim, expanded, flat)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/holoviews/core/data/xarray.py", line 375, in values
    data = data.data
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/core/dataarray.py", line 625, in data
    return self.variable.data
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/core/variable.py", line 347, in data
    return self.values
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/core/variable.py", line 520, in values
    return _as_array_or_item(self._data)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/core/variable.py", line 262, in _as_array_or_item
    data = data.get() if isinstance(data, cupy_array_type) else np.asarray(data)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/numpy/core/_asarray.py", line 102, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/core/indexing.py", line 701, in __array__
    self._ensure_cached()
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/core/indexing.py", line 698, in _ensure_cached
    self.array = NumpyIndexingAdapter(np.asarray(self.array))
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/numpy/core/_asarray.py", line 102, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/core/indexing.py", line 671, in __array__
    return np.asarray(self.array, dtype=dtype)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/numpy/core/_asarray.py", line 102, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/core/indexing.py", line 572, in __array__
    return np.asarray(array[self.key], dtype=None)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/numpy/core/_asarray.py", line 102, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/coding/variables.py", line 70, in __array__
    return self.func(self.array)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/coding/variables.py", line 217, in _scale_offset_decoding
    data = np.array(data, dtype=dtype, copy=True)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/core/indexing.py", line 572, in __array__
    return np.asarray(array[self.key], dtype=None)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/backends/netCDF4_.py", line 91, in __getitem__
    return indexing.explicit_indexing_adapter(
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/core/indexing.py", line 863, in explicit_indexing_adapter
    result = raw_indexing_method(raw_key.tuple)
  File "/home/xin/miniconda3/envs/python38/lib/python3.8/site-packages/xarray/backends/netCDF4_.py", line 114, in _getitem
    raise IndexError(msg)
IndexError: The indexing operation you are attempting to perform is not valid on netCDF4.Variable object. Try loading your data into memory first by calling .load().

The problem with tile sources are that they are kind of unaware of where to start the extents. If you use wheel_zoom + zoom out or set coastline.opts(global_extent=True) or coastline.opts(xlim=(-50, 50), ylim=(-50, 50)), you’ll see something

Also, use active_tools in .opts(), not instantiation.