Interactive plotting with HoloViews

Authors & Contributors

Authors

• Pier Lorenzo Marasco, Provare LTD (UK), @pl-marasco

Contributors

• Alejandro Coca-Castro, The Alan Turing Institute (United Kingdom), @acocac

• Anne Fouilloux, Simula Research Laboratory (Norway), @annefou

Overview

Questions

• How to clip Xarray DataSet with a shapefile to plot an area of interest?
• How to generate interactive maps with Holoviews?
• How to browse multiple times with Holoviews interactive maps?
• How to generate an interactive 1D plot with Holoviews?

Objectives

• Learn about Holoviews and Xarray for creating interactive plots
• Learn to create interactive maps with Holoviews
• Learn to create multiple interactive plots with holoviews
• Learn to create a 1D line plot with holoviews

Context

We will be using HoloViews, a tool part of the HoloViz ecosystem, with Xarray to visualize the Vegetation Condition Index (VCI) [Kog95], a well-established indicator to estimate droughts from remote sensing data.

Data

Data will be downloaded from zenodo using pooch, a very handy python-based library to download and cache your data files locally (see further info here).

import pooch

cgls_file = pooch.retrieve(
    url="https://zenodo.org/record/6969999/files/C_GLS_NDVI_20220101_20220701_Lombardia_S3_2_masked.nc",
    known_hash="md5:be3f16913ebbdb4e7af227f971007b22",
    path=f".",
)

Setup

This episode uses the following main Python packages:

• xarray [HH17] with netCDF4 and h5netcdf engines

• pooch [USR+20]

• rioxarray [SBR+22]

• matplotlib [Hun07]

• cartopy [MetOffice15]

• hvplot [RSB+20]

• geopandas [JdBF+20]

Please install these packages if not already available in your Python environment.

Packages

In this episode, Python packages are imported when we start to use them. However, for best software practices, we recommend you to install and import all the necessary libraries at the top of your Jupyter notebook.

import xarray as xr

Open local dataset

cgls_ds = xr.open_dataset(cgls_file)

Tip

If you get an error with the previous command, check the input file some_hash-C_GLS_NDVI_20220101_20220701_Lombardia_S3_2_masked.netcdf is downloaded locally and it is in the same directory as your Jupyter Notebook.

cgls_ds

<xarray.Dataset> Size: 96MB
Dimensions:  (time: 20, lon: 984, lat: 612)
Coordinates:
  * time     (time) datetime64[ns] 160B 2022-01-01 2022-01-11 ... 2022-07-11
  * lon      (lon) float64 8kB 8.502 8.505 8.508 8.511 ... 11.42 11.42 11.43
  * lat      (lat) float64 5kB 46.5 46.5 46.49 46.49 ... 44.69 44.69 44.68 44.68
Data variables:
    NDVI     (time, lat, lon) float64 96MB ...

xarray.Dataset

• Dimensions:
  • time: 20
  • lon: 984
  • lat: 612
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2022-01-01 ... 2022-07-11

    standard_name :

    t

    long_name :

    t

    axis :

    T

    array(['2022-01-01T00:00:00.000000000', '2022-01-11T00:00:00.000000000',
           '2022-01-21T00:00:00.000000000', '2022-02-01T00:00:00.000000000',
           '2022-02-11T00:00:00.000000000', '2022-02-21T00:00:00.000000000',
           '2022-03-01T00:00:00.000000000', '2022-03-11T00:00:00.000000000',
           '2022-03-21T00:00:00.000000000', '2022-04-01T00:00:00.000000000',
           '2022-04-11T00:00:00.000000000', '2022-04-21T00:00:00.000000000',
           '2022-05-01T00:00:00.000000000', '2022-05-11T00:00:00.000000000',
           '2022-05-21T00:00:00.000000000', '2022-06-01T00:00:00.000000000',
           '2022-06-11T00:00:00.000000000', '2022-06-21T00:00:00.000000000',
           '2022-07-01T00:00:00.000000000', '2022-07-11T00:00:00.000000000'],
          dtype='datetime64[ns]')

  • lon
    (lon)
    float64
    8.502 8.505 8.508 ... 11.42 11.43

    standard_name :

    projection_x_coordinate

    long_name :

    x coordinate

    units :

    degrees_east

    axis :

    X

    array([ 8.501737,  8.504713,  8.507689, ..., 11.42138 , 11.424356, 11.427332])

  • lat
    (lat)
    float64
    46.5 46.5 46.49 ... 44.68 44.68

    standard_name :

    projection_y_coordinate

    long_name :

    y coordinate

    units :

    degrees_north

    axis :

    Y

    array([46.499803, 46.496827, 46.493851, ..., 44.687303, 44.684327, 44.681351])

• Data variables: (1)
  • NDVI
    (time, lat, lon)
    float64
    ...

    [12044160 values with dtype=float64]

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2022-01-01', '2022-01-11', '2022-01-21', '2022-02-01',
                   '2022-02-11', '2022-02-21', '2022-03-01', '2022-03-11',
                   '2022-03-21', '2022-04-01', '2022-04-11', '2022-04-21',
                   '2022-05-01', '2022-05-11', '2022-05-21', '2022-06-01',
                   '2022-06-11', '2022-06-21', '2022-07-01', '2022-07-11'],
                  dtype='datetime64[ns]', name='time', freq=None))

  • lon
    PandasIndex

    PandasIndex(Index([ 8.501737095238099,  8.504713285714303,  8.507689476190507,
             8.51066566666671,  8.513641857142915,  8.516618047619119,
            8.519594238095323,  8.522570428571527,   8.52554661904773,
            8.528522809523935,
           ...
           11.400546619060794, 11.403522809536998, 11.406499000013202,
           11.409475190489406,  11.41245138096561, 11.415427571441814,
           11.418403761918018, 11.421379952394222, 11.424356142870426,
            11.42733233334663],
          dtype='float64', name='lon', length=984))

  • lat
    PandasIndex

    PandasIndex(Index([46.49980333333342, 46.49682714285723, 46.49385095238104,
           46.49087476190485, 46.48789857142866, 46.48492238095247,
           46.48194619047628, 46.47897000000009, 46.47599380952389,
            46.4730176190477,
           ...
           44.70813666666703, 44.70516047619084, 44.70218428571465,
           44.69920809523846, 44.69623190476227, 44.69325571428608,
           44.69027952380989,  44.6873033333337, 44.68432714285751,
           44.68135095238132],
          dtype='float64', name='lat', length=612))

• Attributes: (0)

Clipping data according to a polygon

One of the basic concepts in GIS is to clip data using a vector geometry. Xarray is not directly capable of dealing with vectors but thanks to Rioxarray that can be easily achieved. Rioxarray extends Xarray with most of the features that Rasterio (GDAL) brings.

Read a shapefile with the Area Of Interest (AOI)

import geopandas as gpd

We define the area of interest using the Global Administrative Unit Layers from ASAP - Anomaly Hotspots of Agricultural Production. GeoPandas, a python-based library extending the capabilities of Pandas to deal with geometry and spatial operations, will help to manage geodata.

The full dataset is available at https://agricultural-production-hotspots.ec.europa.eu/files/gaul1_asap_v04.zip but as the version may change and the dataset is quite big, we also created a tiny file containing information about Italy only for the training course

GAUL = gpd.read_file('Italy.geojson')

Data are organized in a tabular structure. For each element an index, data (made of columns) and a geometry are defined.

Geometries are defined through shapely geometry objects with three different basic classes:

• Points and Multi-Points

• Lines and Multi-Lines

• Polygons and Multi-Polygons

GAUL.head(5)

	asap1_id	name1	name1_shr	name0	asap0_id	name0_shr	km2_tot	km2_crop	km2_range	an_crop	an_range	water_lim	geometry
0	1172	Friuli-venezia Giulia	Friuli	Italy	172	Italy	7722	2995	1111	1	1	0	POLYGON ((12.79517 46.646, 12.80945 46.63482, ...
1	1202	Calabria	Calabria	Italy	172	Italy	15203	7152	2127	1	1	1	POLYGON ((16.42261 40.14416, 16.43624 40.12897...
2	1204	Campania	Campania	Italy	172	Italy	13618	7376	2378	1	1	1	POLYGON ((14.16099 41.48058, 14.16653 41.48039...
3	1207	Abruzzi	Abruzzi	Italy	172	Italy	10797	4699	3300	1	1	1	POLYGON ((13.91523 42.89441, 13.92206 42.85505...
4	1211	Emilia-romagna	Emilia-romagna	Italy	172	Italy	22204	14958	771	1	1	1	POLYGON ((12.49151 43.9178, 12.48442 43.90761,...

In the cell below, we subset the polygon geometry in which the name1 field equals to Lombardia.

AOI_name = 'Lombardia'
AOI = GAUL[GAUL.name1 == AOI_name]
AOI_poly = AOI.geometry
AOI_poly

14    POLYGON ((10.23973 46.62177, 10.25084 46.6111,...
Name: geometry, dtype: geometry

In a second step we set a EPSG:4326 Geodetic coordinate reference system to the polygon geometry. To achieve this we need to rely on rioxarray that extends xarray with the rasterio capabilities. The rio accessor is activated through importing rioxarray as has been done at the top.

cgls_ds.rio.write_crs(4326, inplace=True)

<xarray.Dataset> Size: 96MB
Dimensions:      (time: 20, lon: 984, lat: 612)
Coordinates:
  * time         (time) datetime64[ns] 160B 2022-01-01 2022-01-11 ... 2022-07-11
  * lon          (lon) float64 8kB 8.502 8.505 8.508 8.511 ... 11.42 11.42 11.43
  * lat          (lat) float64 5kB 46.5 46.5 46.49 46.49 ... 44.69 44.68 44.68
    spatial_ref  int64 8B 0
Data variables:
    NDVI         (time, lat, lon) float64 96MB ...

xarray.Dataset

• Dimensions:
  • time: 20
  • lon: 984
  • lat: 612
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2022-01-01 ... 2022-07-11

    standard_name :

    t

    long_name :

    t

    axis :

    T

    array(['2022-01-01T00:00:00.000000000', '2022-01-11T00:00:00.000000000',
           '2022-01-21T00:00:00.000000000', '2022-02-01T00:00:00.000000000',
           '2022-02-11T00:00:00.000000000', '2022-02-21T00:00:00.000000000',
           '2022-03-01T00:00:00.000000000', '2022-03-11T00:00:00.000000000',
           '2022-03-21T00:00:00.000000000', '2022-04-01T00:00:00.000000000',
           '2022-04-11T00:00:00.000000000', '2022-04-21T00:00:00.000000000',
           '2022-05-01T00:00:00.000000000', '2022-05-11T00:00:00.000000000',
           '2022-05-21T00:00:00.000000000', '2022-06-01T00:00:00.000000000',
           '2022-06-11T00:00:00.000000000', '2022-06-21T00:00:00.000000000',
           '2022-07-01T00:00:00.000000000', '2022-07-11T00:00:00.000000000'],
          dtype='datetime64[ns]')

  • lon
    (lon)
    float64
    8.502 8.505 8.508 ... 11.42 11.43

    standard_name :

    projection_x_coordinate

    long_name :

    x coordinate

    units :

    degrees_east

    axis :

    X

    array([ 8.501737,  8.504713,  8.507689, ..., 11.42138 , 11.424356, 11.427332])

  • lat
    (lat)
    float64
    46.5 46.5 46.49 ... 44.68 44.68

    standard_name :

    projection_y_coordinate

    long_name :

    y coordinate

    units :

    degrees_north

    axis :

    Y

    array([46.499803, 46.496827, 46.493851, ..., 44.687303, 44.684327, 44.681351])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    array(0)

• Data variables: (1)
  • NDVI
    (time, lat, lon)
    float64
    ...

    [12044160 values with dtype=float64]

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2022-01-01', '2022-01-11', '2022-01-21', '2022-02-01',
                   '2022-02-11', '2022-02-21', '2022-03-01', '2022-03-11',
                   '2022-03-21', '2022-04-01', '2022-04-11', '2022-04-21',
                   '2022-05-01', '2022-05-11', '2022-05-21', '2022-06-01',
                   '2022-06-11', '2022-06-21', '2022-07-01', '2022-07-11'],
                  dtype='datetime64[ns]', name='time', freq=None))

  • lon
    PandasIndex

    PandasIndex(Index([ 8.501737095238099,  8.504713285714303,  8.507689476190507,
             8.51066566666671,  8.513641857142915,  8.516618047619119,
            8.519594238095323,  8.522570428571527,   8.52554661904773,
            8.528522809523935,
           ...
           11.400546619060794, 11.403522809536998, 11.406499000013202,
           11.409475190489406,  11.41245138096561, 11.415427571441814,
           11.418403761918018, 11.421379952394222, 11.424356142870426,
            11.42733233334663],
          dtype='float64', name='lon', length=984))

  • lat
    PandasIndex

    PandasIndex(Index([46.49980333333342, 46.49682714285723, 46.49385095238104,
           46.49087476190485, 46.48789857142866, 46.48492238095247,
           46.48194619047628, 46.47897000000009, 46.47599380952389,
            46.4730176190477,
           ...
           44.70813666666703, 44.70516047619084, 44.70218428571465,
           44.69920809523846, 44.69623190476227, 44.69325571428608,
           44.69027952380989,  44.6873033333337, 44.68432714285751,
           44.68135095238132],
          dtype='float64', name='lat', length=612))

• Attributes: (0)

Once this has been done we can clip the data with the polygon that has been obtained through geopandas at the beginning of the notebook.

NDVI_AOI = cgls_ds.NDVI.rio.clip(AOI_poly, crs=4326)

Visualize with matplotlib

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

fig = plt.figure(1, figsize=[20, 10])

# We're using cartopy and are plotting in PlateCarree projection 
# (see documentation on cartopy)
ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree())
#ax.set_extent([15.5, 27.5, 36, 41], crs=ccrs.PlateCarree()) # lon1 lon2 lat1 lat2
ax.coastlines(resolution='10m')
ax.gridlines(draw_labels=True)

NDVI_AOI.sel(time='2022-06-01').plot(ax=ax, transform=ccrs.PlateCarree(), cmap="RdYlGn")

# One way to customize your title
plt.title("S3 NDVI over Lombardia", fontsize=18)

Text(0.5, 1.0, 'S3 NDVI over Lombardia')

Visualization with HoloViews

import holoviews as hv
import hvplot.xarray
import hvplot.pandas

Plotting data through the HoloViews back-end thanks to the hvplot that acts as high-level plotting API.

NDVI_AOI.isel(time=0).hvplot(cmap="RdYlGn", width=1000, height=1000)

Let’s have a look at the data distribution.

NDVI_AOI[0].hvplot.hist(cmap="RdYlGn",bins=25, width=800, height=700)

Multi-plots using groupby

To be able to visualize interactively all the different available times, we can use groupby time.

NDVI_AOI.hvplot(groupby ='time', cmap="RdYlGn", width=800, height=700)

We can add a histogram to the visualization.

NDVI_AOI.hvplot(groupby='time', cmap='RdYlGn', width=800, height=700 ).hist()

Plot a time series

NDVI_AOI.sel(lat=45.88, lon=8.63, method='nearest').hvplot(ylim=(-0.08, 0.92))

Key Points

• Rioxarray for clipping over an area of interest
• Interactive plot with HoloViews and Xarray

References

[Kog95]

F.N. Kogan. Application of vegetation index and brightness temperature for drought detection. Advances in Space Research, 15(11):91–100, 1995. Natural Hazards: Monitoring and Assessment Using Remote Sensing Technique. URL: https://www.sciencedirect.com/science/article/pii/027311779500079T, doi:https://doi.org/10.1016/0273-1177(95)00079-T.

Packages citation

[Cra23]

Fabio Crameri. Scientific colour maps. 2023. URL: https://zenodo.org/record/1243862, doi:10.5281/ZENODO.1243862.

[HH17]

S. Hoyer and J. Hamman. Xarray: N-D labeled arrays and datasets in Python. Journal of Open Research Software, 2017. URL: https://doi.org/10.5334/jors.148, doi:10.5334/jors.148.

[Hun07]

J. D. Hunter. Matplotlib: a 2d graphics environment. Computing in Science & Engineering, 9(3):90–95, 2007. doi:10.1109/MCSE.2007.55.

[JdBF+20]

Kelsey Jordahl, Joris Van den Bossche, Martin Fleischmann, Jacob Wasserman, James McBride, Jeffrey Gerard, Jeff Tratner, Matthew Perry, Adrian Garcia Badaracco, Carson Farmer, Geir Arne Hjelle, Alan D. Snow, Micah Cochran, Sean Gillies, Lucas Culbertson, Matt Bartos, Nick Eubank, maxalbert, Aleksey Bilogur, Sergio Rey, Christopher Ren, Dani Arribas-Bel, Leah Wasser, Levi John Wolf, Martin Journois, Joshua Wilson, Adam Greenhall, Chris Holdgraf, Filipe, and François Leblanc. Geopandas/geopandas: v0.8.1. July 2020. URL: https://doi.org/10.5281/zenodo.3946761, doi:10.5281/zenodo.3946761.

[RSB+20]

Philipp Rudiger, Jean-Luc Stevens, James A. Bednar, Bas Nijholt, Andrew, Chris B, Achim Randelhoff, Jon Mease, Vasco Tenner, maxalbert, Markus Kaiser, ea42gh, Jordan Samuels, stonebig, Florian LB, Andrew Tolmie, Daniel Stephan, Scott Lowe, John Bampton, henriqueribeiro, Irv Lustig, Julia Signell, Justin Bois, Leopold Talirz, Lukas Barth, Maxime Liquet, Ram Rachum, Yuval Langer, arabidopsis, and kbowen. Holoviz/holoviews: version 1.13.3. June 2020. URL: https://doi.org/10.5281/zenodo.3904606, doi:10.5281/zenodo.3904606.

[SBR+22]

Alan D. Snow, David Brochart, Martin Raspaud, Ray Bell, RichardScottOZ, Taher Chegini, Alessandro Amici, Rémi Braun, Andrew Annex, David Hoese, Fred Bunt, GBallesteros, Joe Hamman, Markus Zehner, Mauricio Cordeiro, Scott Henderson, Seth Miller, The Gitter Badger, Tom Augspurger, apiwat-chantawibul, and pmallas. Corteva/rioxarray: 0.11.1 release. April 2022. URL: https://doi.org/10.5281/zenodo.6478182, doi:10.5281/zenodo.6478182.

[USR+20]

Leonardo Uieda, Santiago Rubén Soler, Rémi Rampin, Hugo van Kemenade, Matthew Turk, Daniel Shapero, Anderson Banihirwe, and John Leeman. Pooch: a friend to fetch your data files. Journal of Open Source Software, 5(45):1943, 2020. URL: https://doi.org/10.21105/joss.01943, doi:10.21105/joss.01943.

[MetOffice15]

Met Office. Cartopy: a cartographic python library with a matplotlib interface. Exeter, Devon, 2010 - 2015. URL: http://scitools.org.uk/cartopy.