Getting started¶

The xdatasets library enables users to effortlessly access a vast collection of earth observation datasets that are compatible with xarray formats.

The library adopts an opinionated approach to data querying and caters to the specific needs of certain user groups, such as hydrologists, climate scientists, and engineers. One of the functionalities of xdatasets is the ability to extract data at a specific location or within a designated region, such as a watershed or municipality, while also enabling spatial and temporal operations.

To use xdatasets, users must employ a query. For instance, a straightforward query to extract the variables t2m (2m temperature) and tp (Total precipitation) from the era5_reanalysis_single_levels dataset at two geographical positions (Montreal and Toronto) could be as follows:

query = {
    "datasets": {"era5_reanalysis_single_levels": {'variables': ["t2m", "tp"]}},
    "space": {
        "clip": "point", # bbox, point or polygon
        "geometry": {'Montreal' : (45.508888, -73.561668),
                     'Toronto' : (43.651070, -79.347015)
                    }
    }
}

An example of a more complex query would look like the one below.

Note Don’t worry! Below, you’ll find additional examples that will assist in understanding each parameter in the query, as well as the possible combinations.

This query calls the same variables as above. However, instead of specifying geographical positions, a GeoPandas.DataFrame is used to provide features (such as shapefiles or geojson) for extracting data within each of them. Each polygon is identified using the unique identifier Station, and a spatial average is computed within each one (aggregation: True). The dataset, initially at an hourly time step, is converted into a daily time step while applying one or more temporal aggregations for each variable as prescribed in the query. xdatasets ultimately returns the dataset for the specified date range and time zone.

query = {
    "datasets": {"era5_reanalysis_single_levels": {'variables': ["t2m", "tp"]}},
    "space": {
        "clip": "polygon", # bbox, point or polygon
        "averaging": True, # spatial average of the variables within each polygon
        "geometry": gdf,
        "unique_id": "Station" # unique column name in geodataframe
    },
    "time": {
        "timestep": "D",
        "aggregation": {"tp": np.nansum,
                        "t2m": [np.nanmax, np.nanmin]},

        "start": '2000-01-01',
        "end": '2020-05-31',
        "timezone": 'America/Montreal',
    },
}

Query climate datasets¶

In order to use xdatasets, you must import at least xdatasets, pandas, geopandas, and numpy. Additionally, we import pathlib to interact with files.

[1]:

import os
import warnings
from pathlib import Path

warnings.simplefilter("ignore")

os.environ["USE_PYGEOS"] = "0"
import geopandas as gpd

# Visualization
import hvplot.pandas  # noqa
import hvplot.xarray  # noqa-
import numpy as np
import pandas as pd
import panel as pn  # noqa

import xdatasets as xd

Clip on polygons with no averaging in space¶

First, let’s explore specific polygon features. With xdatasets, you can access geographical datasets, such as watershed boundaries linked to streamflow stations. These datasets follow a nomenclature where they are named after the hydrological dataset, with "_polygons" appended. For example, if the hydrological dataset is named deh, its corresponding watershed boundaries dataset will be labeled deh_polygons. The query below retrieves all polygons for the deh_polygons dataset.

gdf = xd.Query(
    **{
        "datasets": "deh_polygons"
}).data

gdf

As the data is loaded into memory, the process of loading all polygons may take some time. To expedite this, we recommend employing filters, as illustrated below. It’s important to note that the filters are consistent for both hydrological and corresponding geographical datasets. Consequently, only watershed boundaries associated with existing hydrological data will be returned.

[7]:

import xdatasets as xd

gdf = xd.Query(
    **{
        "datasets": {
            "deh_polygons": {
                "id": ["0421*"],
            }
        }
    }
).data.reset_index()

gdf

[7]:

	Station	Superficie	geometry
0	042102	623.479187	POLYGON ((-78.5712 46.70742, -78.57112 46.7073...
1	042103	579.479614	POLYGON ((-78.49014 46.64514, -78.4901 46.6450...

Let’s examine the geographic locations of the polygon features.

[8]:

gdf.hvplot(
    geo=True,
    tiles="ESRI",
    color="Station",
    alpha=0.8,
    width=750,
    height=450,
    legend="top",
    hover_cols=["Station", "Superficie"],
)

---------------------------------------------------------------------------
DataError                                 Traceback (most recent call last)
Cell In[8], line 1
----> 1 gdf.hvplot(
      2     geo=True,
      3     tiles="ESRI",
      4     color="Station",
      5     alpha=0.8,
      6     width=750,
      7     height=450,
      8     legend="top",
      9     hover_cols=["Station", "Superficie"],
     10 )

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/hvplot/plotting/core.py:95, in hvPlotBase.__call__(self, x, y, kind, **kwds)
     92         plot = self._get_converter(x, y, kind, **kwds)(kind, x, y)
     93         return pn.panel(plot, **panel_dict)
---> 95 return self._get_converter(x, y, kind, **kwds)(kind, x, y)

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/hvplot/converter.py:1645, in HoloViewsConverter.__call__(self, kind, x, y)
   1643     columns = [c for c in data.columns if c != 'geometry']
   1644     shape_dims = ['Longitude', 'Latitude'] if self.geo else ['x', 'y']
-> 1645     dataset = Dataset(data, kdims=shape_dims + columns)
   1646 elif self.datatype == 'xarray':
   1647     import xarray as xr

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/holoviews/core/data/__init__.py:333, in Dataset.__init__(self, data, kdims, vdims, **kwargs)
    331 (data, self.interface, dims, extra_kws) = initialized
    332 super().__init__(data, **dict(kwargs, **dict(dims, **extra_kws)))
--> 333 self.interface.validate(self, validate_vdims)
    335 # Handle _pipeline property
    336 if input_pipeline is None:

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/holoviews/core/data/pandas.py:166, in PandasInterface.validate(cls, dataset, vdims)
    164 not_found = [d for d in dimensions if d not in cols]
    165 if not_found:
--> 166     raise DataError("Supplied data does not contain specified "
    167                     "dimensions, the following dimensions were "
    168                     "not found: %s" % repr(not_found), cls)

DataError: Supplied data does not contain specified dimensions, the following dimensions were not found: ['Longitude', 'Latitude']

PandasInterface expects tabular data, for more information on supported datatypes see http://holoviews.org/user_guide/Tabular_Datasets.html

The following query seeks the variables t2m and tp from the era5_reanalysis_single_levels dataset, covering the period between January 1, 1959, and September 30, 1961, for the three polygons mentioned earlier. It is important to note that as aggregation is set to False, no spatial averaging will be conducted, and a mask (raster) will be returned for each polygon.

[9]:

query = {
    "datasets": {"era5_reanalysis_single_levels": {"variables": ["t2m", "tp"]}},
    "space": {
        "clip": "polygon",  # bbox, point or polygon
        "averaging": False,  # spatial average of the variables within each polygon
        "geometry": gdf,
        "unique_id": "Station",  # unique column name in geodataframe
    },
    "time": {
        "start": "1959-01-01",
        "end": "1961-08-31",
    },
}

xds = xd.Query(**query)

By accessing the data attribute, you can view the data obtained from the query. For each variable, the dimensions of time, latitude, longitude, and Station (the unique ID) are included. In addition, there is another variable called weights that is returned. This variable specifies the weight that should be assigned to each pixel if spatial averaging is conducted over a mask (polygon).

[10]:

xds.data

Weights are much easier to comprehend visually, so let’s examine the weights returned for the station 042102. Notice that when selecting a single feature (Station 042102 in this case), the shape of our spatial dimensions is reduced to a 3x2 pixel area (longitude x latitude) that encompasses the entire feature.

[11]:

station = "042102"

ds_station = xds.data.sel(Station=station)
ds_clipped = xds.bbox_clip(ds_station).squeeze()
ds_clipped

[12]:

(
    (
        ds_clipped.t2m.isel(time=0).hvplot(
            title="The 2m temperature for pixels that intersect with the polygon on January 1, 1959",
            tiles="ESRI",
            geo=True,
            alpha=0.6,
            colormap="isolum",
            width=750,
            height=450,
        )
        * gdf[gdf.Station == station].hvplot(
            geo=True,
            width=750,
            height=450,
            legend="top",
            hover_cols=["Station", "Superficie"],
        )
    )
    + ds_clipped.weights.hvplot(
        title="The weights that should be assigned to each pixel when performing spatial averaging",
        tiles="ESRI",
        alpha=0.6,
        colormap="isolum",
        geo=True,
        width=750,
        height=450,
    )
    * gdf[gdf.Station == station].hvplot(
        geo=True,
        width=750,
        height=450,
        legend="top",
        hover_cols=["Station", "Superficie"],
    )
).cols(1)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[12], line 3
      1 (
      2     (
----> 3         ds_clipped.t2m.isel(time=0).hvplot(
      4             title="The 2m temperature for pixels that intersect with the polygon on January 1, 1959",
      5             tiles="ESRI",
      6             geo=True,
      7             alpha=0.6,
      8             colormap="isolum",
      9             width=750,
     10             height=450,
     11         )
     12         * gdf[gdf.Station == station].hvplot(
     13             geo=True,
     14             width=750,
     15             height=450,
     16             legend="top",
     17             hover_cols=["Station", "Superficie"],
     18         )
     19     )
     20     + ds_clipped.weights.hvplot(
     21         title="The weights that should be assigned to each pixel when performing spatial averaging",
     22         tiles="ESRI",
     23         alpha=0.6,
     24         colormap="isolum",
     25         geo=True,
     26         width=750,
     27         height=450,
     28     )
     29     * gdf[gdf.Station == station].hvplot(
     30         geo=True,
     31         width=750,
     32         height=450,
     33         legend="top",
     34         hover_cols=["Station", "Superficie"],
     35     )
     36 ).cols(1)

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/hvplot/plotting/core.py:95, in hvPlotBase.__call__(self, x, y, kind, **kwds)
     92         plot = self._get_converter(x, y, kind, **kwds)(kind, x, y)
     93         return pn.panel(plot, **panel_dict)
---> 95 return self._get_converter(x, y, kind, **kwds)(kind, x, y)

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/hvplot/converter.py:1671, in HoloViewsConverter.__call__(self, kind, x, y)
   1668     obj = gv.project(obj, projection=self.output_projection)
   1670 if not (self.datashade or self.rasterize or self.downsample):
-> 1671     layers = self._apply_layers(obj)
   1672     layers = _transfer_opts_cur_backend(layers)
   1673     return layers

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/hvplot/converter.py:1826, in HoloViewsConverter._apply_layers(self, obj)
   1824     tiles = self._get_tiles(self.tiles)
   1825 else:
-> 1826     tiles = self._get_tiles(self.tiles, lib='geoviews')
   1827 tiles = tiles.opts(clone=True, **self.tiles_opts)
   1828 obj = tiles * obj

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/hvplot/converter.py:1868, in HoloViewsConverter._get_tiles(self, source, lib)
   1858 if tiles is None:
   1859     msg = (
   1860         f'{tile_source} tiles not recognized. tiles must be either True, a '
   1861         'xyzservices.TileProvider instance, a HoloViews'
   (...)   1866         + '.'
   1867     )
-> 1868     raise ValueError(msg)
   1869 return tiles

ValueError: EsriImagery tiles not recognized. tiles must be either True, a xyzservices.TileProvider instance, a HoloViews or Geoviews basemap string (one of CartoDark, CartoEco, CartoLight, CartoMidnight, ESRI, OSM, StamenLabels, StamenTerrain, StamenTerrainRetina, StamenToner, StamenTonerBackground, StamenWatercolor, Wikipedia), a HoloViews Tiles instance, a Geoviews WMTS instance.

The two plots depicted above show the 2m temperature for each pixel that intersects with the polygon from Station 042102 and the corresponding weights to be applied to each pixel. In the lower plot, it is apparent that the majority of the polygon is situated in the central pixels, which results in those pixels having a weight of approximately 80%. It is evident that the two lower and the upper pixels have much less intersection with the polygon, which results in their respective weights being smaller (hover on the plot to verify the weights).

In various libraries, either all pixels that intersect with the geometries are kept, or only pixels with centers within the polygon are retained. However, as shown in the previous example, utilizing such methods can introduce significant biases in the final calculations.

Query hydrological datasets¶

Hydrological queries are still being tested and output format is likely to change. Stay tuned!

[21]:

query = {"datasets": "deh"}
xds = xd.Query(**query)
xds.data

[21]:

<xarray.Dataset> Size: 1GB
Dimensions:        (id: 745, time: 60631, variable: 2, spatial_agg: 2,
                    timestep: 1, time_agg: 1, source: 1)
Coordinates: (12/15)
  * id             (id) object 6kB '010101' '010801' ... '104804' '120201'
  * time           (time) datetime64[ns] 485kB 1860-01-01 ... 2025-12-31
  * variable       (variable) object 16B 'level' 'streamflow'
  * spatial_agg    (spatial_agg) object 16B 'point' 'watershed'
  * timestep       (timestep) object 8B 'D'
  * time_agg       (time_agg) object 8B 'mean'
    ...             ...
    latitude       (id) float32 3kB 48.48 48.1 48.19 48.19 ... 56.78 55.86 46.97
    longitude      (id) float32 3kB -64.53 -65.46 -65.56 ... -64.88 -70.86
    name           (id) object 6kB 'La Grande Rivière' ... 'Dauphine'
    province       (id) object 6kB 'QC' 'QC' 'QC' 'QC' ... 'QC' 'QC' 'QC' 'QC'
    regulated      (id) object 6kB 'Natural' 'Natural' ... 'Natural' 'Natural'
    start_date     (variable, id, spatial_agg, timestep, time_agg, source) datetime64[ns] 24kB ...
Data variables:
    level          (id, time, variable, spatial_agg, timestep, time_agg, source) float32 723MB dask.array<chunksize=(1, 60631, 1, 1, 1, 1, 1), meta=np.ndarray>
    streamflow     (id, time, variable, spatial_agg, timestep, time_agg, source) float32 723MB dask.array<chunksize=(1, 60631, 1, 1, 1, 1, 1), meta=np.ndarray>

[22]:

ds = (
    xd.Query(
        **{
            "datasets": {
                "deh": {
                    "id": ["020*"],
                    "regulated": ["Natural"],
                    "variables": ["streamflow"],
                }
            },
            "time": {"start": "1970-01-01", "minimum_duration": (10 * 365, "d")},
        }
    )
    .data.squeeze()
    .load()
)

ds

[23]:

query = {"datasets": "hydat"}
xds = xd.Query(**query)
xds.data

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[23], line 2
      1 query = {"datasets": "hydat"}
----> 2 xds = xd.Query(**query)
      3 xds.data

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/xdatasets/core.py:130, in Query.__init__(self, datasets, space, time, catalog_path)
    127 self.space = self._resolve_space_params(**space)
    128 self.time = self._resolve_time_params(**time)
--> 130 self.load_query(datasets=self.datasets, space=self.space, time=self.time)

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/xdatasets/core.py:258, in Query.load_query(self, datasets, space, time)
    255     except:  # noqa: S110
    256         pass
--> 258     ds_one = self._process_one_dataset(
    259         dataset_name=dataset_name,
    260         variables=variables_name,
    261         space=space,
    262         time=time,
    263         **kwargs,
    264     )
    265     dsets.append(ds_one)
    267 try:
    268     # Try naively merging datasets into single dataset

File ~/checkouts/readthedocs.org/user_builds/xdatasets/conda/stable/lib/python3.12/site-packages/xdatasets/core.py:301, in Query._process_one_dataset(self, dataset_name, variables, space, time, **kwargs)
    298     dataset_category = "user-provided"
    300 elif isinstance(dataset_name, str):
--> 301     dataset_category = [
    302         category for category in self.catalog._entries.keys() for name in self.catalog[category]._entries.keys() if name == dataset_name
    303     ][0]
    305 if dataset_category in ["atmosphere"]:
    306     with warnings.catch_warnings():

IndexError: list index out of range

Getting started¶

Query climate datasets¶

Clip by points (sites)¶

Clip on polygons with no averaging in space¶

Clip on polygons with averaging in space¶

Bounding box (bbox) around polygons¶

Query hydrological datasets¶