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
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 ; 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
``````

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)
``````

1 Like

@ahuang11 Wow, thanks! That’s quite nice! I will play with my data using different tiles
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:

``````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
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__