Contents

Interactive maps

In this notebook, we will explore other capabilities of the hvPlot ecosystem. The associated functions allow to create plots with many interactive tools and features, including linked panning and brushing, and hover inspectors.

In particular, we will look at the VectorField method, using the AIMS eReefs database.

See also

For other inspirational functionalities, you might be interested in the examples from the HoloViews and GeoViews gallery.

Load the required Python libraries

First of all, load the necessary libraries. These are the ones we discussed previously:

  • numpy

  • matplotlib

  • cartopy

  • xarray

  • holoviews

  • geoviews

import os
import numpy as np
import xarray as xr

import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
cartopy.config['data_dir'] = os.getenv('CARTOPY_DIR', cartopy.config.get('data_dir'))

import cmocean

import holoviews as hv
from holoviews import opts, dim

import geoviews as gv
import geoviews.feature as gf
from cartopy import crs

gv.extension('bokeh')

import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)

Build a multi-file dataset

We will use the open_mfdataset function from xArray to open multiple netCDF files into a single xarray Dataset.

We will query load the GBR4km dataset from the AIMS server, so let’s first define the base URL.

base_url = "http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-"

For the sake of the demonstration, we will only use 1 month:

month_st = 1   # Starting month 
month_ed = 1   # Ending month 
year = 2018    # Year

# Based on the server the file naming convention 
hydrofiles = [f"{base_url}{year}-{month:02}.nc" for month in range(month_st, month_ed+1)]
hydrofiles

['http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-2018-01.nc']

Loading the dataset into xArray

Using xArray, we open these files into a Dataset:

ds_hydro = xr.open_mfdataset(hydrofiles)
ds_hydro

<xarray.Dataset>
Dimensions:      (k: 17, latitude: 723, longitude: 491, time: 31)
Coordinates:
    zc           (k) float64 dask.array<chunksize=(17,), meta=np.ndarray>
  * time         (time) datetime64[ns] 2017-12-31T14:00:00 ... 2018-01-30T14:...
  * latitude     (latitude) float64 -28.7 -28.67 -28.64 ... -7.096 -7.066 -7.036
  * longitude    (longitude) float64 142.2 142.2 142.2 ... 156.8 156.8 156.9
Dimensions without coordinates: k
Data variables:
    mean_cur     (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    salt         (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    temp         (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    u            (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    v            (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    mean_wspeed  (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    eta          (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    wspeed_u     (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    wspeed_v     (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
Attributes: (12/19)
    Conventions:                     CF-1.0
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T12:50:07+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    aims_ncaggregate_firstDate:      2018-01-01T00:00:00+10:00
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time

Clip the Dataset

To reduce the Dataset size we will clip the spatial extent based on longitudinal and latitudinal values.

As we already saw this is easely done using the sel function with the slice method.

min_lon = 146     # lower left longitude
min_lat = -21     # lower left latitude
max_lon = 150     # upper right longitude
max_lat = -16     # upper right latitude

# Defining the boundaries
lon_bnds = [min_lon, max_lon]
lat_bnds = [min_lat, max_lat]


# Performing the reduction 
ds_hydro_clip = ds_hydro.sel(latitude=slice(*lat_bnds), longitude=slice(*lon_bnds))
ds_hydro_clip

<xarray.Dataset>
Dimensions:      (k: 17, latitude: 167, longitude: 134, time: 31)
Coordinates:
    zc           (k) float64 dask.array<chunksize=(17,), meta=np.ndarray>
  * time         (time) datetime64[ns] 2017-12-31T14:00:00 ... 2018-01-30T14:...
  * latitude     (latitude) float64 -20.99 -20.96 -20.93 ... -16.04 -16.01
  * longitude    (longitude) float64 146.0 146.0 146.1 ... 149.9 150.0 150.0
Dimensions without coordinates: k
Data variables:
    mean_cur     (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 167, 134), meta=np.ndarray>
    salt         (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 167, 134), meta=np.ndarray>
    temp         (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 167, 134), meta=np.ndarray>
    u            (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 167, 134), meta=np.ndarray>
    v            (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 167, 134), meta=np.ndarray>
    mean_wspeed  (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    eta          (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    wspeed_u     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    wspeed_v     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
Attributes: (12/19)
    Conventions:                     CF-1.0
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T12:50:07+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    aims_ncaggregate_firstDate:      2018-01-01T00:00:00+10:00
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time

We will use for this example the surface values (i.e. the last z-coordinate at -0.5 m), this is done by using isel function on the Xarray Dataset:

zcvar = -1  # last z-coordinate at 0.5 m

new_ds = ds_hydro_clip.isel(k=zcvar)
new_ds

<xarray.Dataset>
Dimensions:      (latitude: 167, longitude: 134, time: 31)
Coordinates:
    zc           float64 dask.array<chunksize=(), meta=np.ndarray>
  * time         (time) datetime64[ns] 2017-12-31T14:00:00 ... 2018-01-30T14:...
  * latitude     (latitude) float64 -20.99 -20.96 -20.93 ... -16.04 -16.01
  * longitude    (longitude) float64 146.0 146.0 146.1 ... 149.9 150.0 150.0
Data variables:
    mean_cur     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    salt         (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    temp         (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    u            (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    v            (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    mean_wspeed  (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    eta          (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    wspeed_u     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    wspeed_v     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
Attributes: (12/19)
    Conventions:                     CF-1.0
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T12:50:07+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    aims_ncaggregate_firstDate:      2018-01-01T00:00:00+10:00
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time

Groupby function

Using groupby with time coordinates, we take the maximum of each Data variables over a 2 weeks timeframe (‘2W’). We also reduce the Dataset to the following variables:

  • ‘mean_cur’,

  • ‘u’, and

  • ‘v’

week_ds = new_ds.resample(time='2W').max(dim='time').drop('zc')[['mean_cur','u','v']]
week_ds

<xarray.Dataset>
Dimensions:    (latitude: 167, longitude: 134, time: 4)
Coordinates:
  * time       (time) datetime64[ns] 2017-12-31 2018-01-14 2018-01-28 2018-02-11
  * latitude   (latitude) float64 -20.99 -20.96 -20.93 ... -16.07 -16.04 -16.01
  * longitude  (longitude) float64 146.0 146.0 146.1 146.1 ... 149.9 150.0 150.0
Data variables:
    mean_cur   (time, latitude, longitude) float32 dask.array<chunksize=(1, 167, 134), meta=np.ndarray>
    u          (time, latitude, longitude) float32 dask.array<chunksize=(1, 167, 134), meta=np.ndarray>
    v          (time, latitude, longitude) float32 dask.array<chunksize=(1, 167, 134), meta=np.ndarray>

We then load the newly created Xarray Dataset in memory:

load_week = week_ds.load()

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-8-e657c306f6eb> in <module>
----> 1 load_week = week_ds.load()

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/xarray/core/dataset.py in load(self, **kwargs)
    796 
    797             # evaluate all the dask arrays simultaneously
--> 798             evaluated_data = da.compute(*lazy_data.values(), **kwargs)
    799 
    800             for k, data in zip(lazy_data, evaluated_data):

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/dask/base.py in compute(*args, **kwargs)
    563         postcomputes.append(x.__dask_postcompute__())
    564 
--> 565     results = schedule(dsk, keys, **kwargs)
    566     return repack([f(r, *a) for r, (f, a) in zip(results, postcomputes)])
    567 

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/dask/threaded.py in get(dsk, result, cache, num_workers, pool, **kwargs)
     74                 pools[thread][num_workers] = pool
     75 
---> 76     results = get_async(
     77         pool.apply_async,
     78         len(pool._pool),

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/dask/local.py in get_async(apply_async, num_workers, dsk, result, cache, get_id, rerun_exceptions_locally, pack_exception, raise_exception, callbacks, dumps, loads, **kwargs)
    474             # Main loop, wait on tasks to finish, insert new ones
    475             while state["waiting"] or state["ready"] or state["running"]:
--> 476                 key, res_info, failed = queue_get(queue)
    477                 if failed:
    478                     exc, tb = loads(res_info)

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/dask/local.py in queue_get(q)
    131 
    132     def queue_get(q):
--> 133         return q.get()
    134 
    135 

/usr/share/miniconda/envs/envireef/lib/python3.8/queue.py in get(self, block, timeout)
    168             elif timeout is None:
    169                 while not self._qsize():
--> 170                     self.not_empty.wait()
    171             elif timeout < 0:
    172                 raise ValueError("'timeout' must be a non-negative number")

/usr/share/miniconda/envs/envireef/lib/python3.8/threading.py in wait(self, timeout)
    300         try:    # restore state no matter what (e.g., KeyboardInterrupt)
    301             if timeout is None:
--> 302                 waiter.acquire()
    303                 gotit = True
    304             else:

KeyboardInterrupt:

Magnitude and angle

HvPlot vectorfield accepts 2d arrays of magnitude and angle on a grid and produces an array of vectors.

See also

A GeoViews example of VectorField declaration.

We will create these variables by taking advantage of Xarray DataArrays and Datasets arithmetic operators and numpy array functions:

# Convert current u and v to magnitude and angle
mag = np.sqrt(load_week.u**2 + load_week.v**2)
angle = (np.pi/2.) - np.arctan2(load_week.u/mag, load_week.v/mag)

We now add these 2 variables to our Xarray Datasets:

load_week["mag"] = (['time', 'latitude', 'longitude'],  mag)
load_week["angle"] = (['time', 'latitude', 'longitude'],  angle)

load_week

<xarray.Dataset>
Dimensions:    (latitude: 167, longitude: 134, time: 4)
Coordinates:
  * time       (time) datetime64[ns] 2017-12-31 2018-01-14 2018-01-28 2018-02-11
  * latitude   (latitude) float64 -20.99 -20.96 -20.93 ... -16.07 -16.04 -16.01
  * longitude  (longitude) float64 146.0 146.0 146.1 146.1 ... 149.9 150.0 150.0
Data variables:
    mean_cur   (time, latitude, longitude) float32 nan nan nan ... 0.2468 0.2474
    u          (time, latitude, longitude) float32 nan nan ... -0.08577 -0.07427
    v          (time, latitude, longitude) float32 nan nan ... 0.09289 0.08283
    mag        (time, latitude, longitude) float32 nan nan nan ... 0.1264 0.1113
    angle      (time, latitude, longitude) float32 nan nan nan ... 2.316 2.302

Dataset temporal visualisation

Using GeoViews we visualise the weekly max surface magnitude:

# Specify the dataset, its coordinates and requested variable 
dataset = gv.Dataset(load_week, ['longitude', 'latitude', 'time'], 
                     'mag', crs=crs.PlateCarree())
images = dataset.to(gv.Image,dynamic=True)

# Loading coastlines from Cartopy
coastline = gf.coastline(line_width=2,line_color='k').opts(projection=ccrs.PlateCarree(),scale='10m')

# Loading land mask from Cartopy
land = gf.land.options(scale='10m', fill_color='lightgray')

# Slider location
hv.output(widget_location='bottom')

# Create stack of images grouped by time
images.opts(active_tools=['wheel_zoom', 'pan'], cmap=cmocean.cm.speed,
            colorbar=True, width=500, height=500, clim=(0.1,1.0)) * coastline * land

Adding a vector field

To plot the vector field we will first resample the dataset and take one point every 7 times. It allows to have less arrows on the map making it more readable…

resample = load_week.isel(longitude=slice(None, None, 7), latitude=slice(None, None, 7))
resample

<xarray.Dataset>
Dimensions:    (latitude: 24, longitude: 20, time: 4)
Coordinates:
  * time       (time) datetime64[ns] 2017-12-31 2018-01-14 2018-01-28 2018-02-11
  * latitude   (latitude) float64 -20.99 -20.78 -20.57 ... -16.58 -16.37 -16.16
  * longitude  (longitude) float64 146.0 146.2 146.4 146.6 ... 149.6 149.8 150.0
Data variables:
    mean_cur   (time, latitude, longitude) float32 nan nan nan ... 0.2565 0.142
    u          (time, latitude, longitude) float32 nan nan ... -0.05966 -0.06938
    v          (time, latitude, longitude) float32 nan nan ... -0.009086 0.04427
    mag        (time, latitude, longitude) float32 nan nan ... 0.06035 0.0823
    angle      (time, latitude, longitude) float32 nan nan nan ... 3.293 2.574
# Maximum magnitude value
resample.mag.max().item()

0.9817500114440918

# x and y 1d coordinates
lat = resample.latitude
lon = resample.longitude

# Number of time step (4 based on the groupby approach 2W)
nb = resample.angle.shape[0]

# We will normalise the arrow to avoid changes in scale as the
# time evolve...
max_mag = resample.mag.max()

# Create a disctionary of VectorField values at each time interval
a_var = resample.time.values
y_list = []
for i in range(nb):
    y_list.append(gv.VectorField((lon, lat, resample.angle[i,:,:],
                                  resample.mag[i,:,:]/max_mag), 
                                 crs=crs.PlateCarree()))
dict_y = {a_var[i]:y_list[i] for i in range(nb)}

# create HoloMap object 
hmap = hv.HoloMap(dict_y, kdims="time").opts(opts.VectorField(magnitude=dim('Magnitude')*0.25, 
                                                              color='k', 
                                                              width=600, height=500,
                                                              pivot='tip', line_width=1, 
                                                              title='eReefs GBR 4km', 
                                                              rescale_lengths=False,
                                                              projection=crs.PlateCarree()))

hmap

Ok let’s combine everything together:

images.opts(active_tools=['wheel_zoom', 'pan'], cmap=cmocean.cm.speed,
            colorbar=True, width=500, height=500, clim=(0.1,1.0)) * coastline * land * hmap

The image can be saved as a html file by uncomenting the cell below:

# figs = images.opts(cmap=cmocean.cm.speed, colorbar=True, width=500, height=500) * coastline * land * hmap
# hv.save(figs, 'graph.html')

See also

To see the above example in full interactivity, you might want to use the following html link

# GOTO: https://raw.githack.com
# https://github.com/TristanSalles/EnviReef/main/book/5-xarray/examples/graph.html
# from IPython.display import IFrame
# from IPython.core.display import display
# display(IFrame('https://raw.githack.com/TristanSalles/EnviReef/main/book/5-xarray/examples/graph.html', width='950px', height='500px'))
Xarray maps Box streaming