Hyrax 1.18.0 is now available.
Read the release notes →

  •  Icon
  •  Icon
  •  Icon
  •  Icon

What OPeNDAP data is needed to create an Xarray Dataset? Multiple URLs case

What OPeNDAP data is needed to create an Xarray Dataset? Multiple URLs case
Creating an aggregation can rapidly become an obstacle for scientific discovery due to some safety features. Why? Here we describe what is going on in the background.

Creating an aggregate Xarray Dataset from N remote OPeNDAP URLs can turn a simple first step into a serious bottleneck for scientific discovery. That first step can often be surprisingly slow, slow enough to become a real barrier to data access. Why? Several culprits are usually in play: a data layout that isn’t suited to aggregation, authentication issues, a flaky internet connection, and so on. Experienced Xarray users tend to spot these quickly. But there’s a subtler cause that, in my experience, is actually the most common culprit behind slow Xarray+OPeNDAP workflows but it is the hardest to diagnose: Xarray’s own built-in caution. Left to its defaults, it issues far too many GET requests to build the aggregation — NOT a bug, but a safety behavior that a user who understands what’s happening underneath can avoid.

Even when the remote OPeNDAP server scales to meet the demands, producing too many GET requests during aggregation can become unnecessarily slow. Below we demonstrate a very simple case that can trigger too many (unnecessary) requests, even with the simplest data layout: a Level 3 dataset.

How many requests?

The performance a user experiences when opening an Xarray dataset is directly related to the answer to this question. So, how many requests are made to the remote OPeNDAP server? When aggregating N remote datasets, the following requests are fetched, each with its unique OPeNDAP URL and, importantly, latency:

  • N metadata requests (i.e. the DMR, each containing all the metadata of each file, such as variable names, attributes, description, etc)
  • N dimension data requests (each containing all the 1D array dimensions from each remote file).
  • M*N requests, one per non-dimension coordinate array from each URL. An example of a non-dimension coordinate array is a 2D Latitude array.

As described above, given N remote URLs, there can be at least N*(2+M) total requests done to the server, possibly more, when aggregating via Xarray. Since there is no restriction on how big a coordinate array (such as a 2D Latitude) can be, then “opening” an Xarray Dataset can be so slow that it effectively becomes a barrier for scientific analysis. In the case of Level 3 data (where the coordinates are identical across all remote files), such amount of downloads is completely unnecessary.

Below we demonstrate this with a very simple Level 3 dataset from DAYMET. On a separate blog post we will provide a workaround making use of the DAP4 protocol.

Tutorial

This tutorial will demonstrate:

  1. How to install the necessary libraries.
  2. How to discover OPeNDAP urls from NASA’s Earthdata Cloud.
  3. How to Authenticate in order to stream data.
  4. How to create an Xarray Dataset.
  5. How to verify the requests made to the remote OPeNDAP server.

Create an isolated environment

Here, we open a terminal and create an isolated containerized environment, to ensure all the work below is reproducible. We will use Python 3.12, and available libraries for such environment. And we will use a jupyterlab environment as our working station.

Terminal
$ mamba create -n opendap_env -c conda-forge python=3.12 ipython pydap jupyterlab earthaccess netCDF4 dask
$ mamba activate opendap_env
$ jupyter lab

Assuming the jupyter notebook session is up and running, we begin by importing the necessary libraries and methods

Python
import xarray as xr
import datetime as dt
import earthaccess

# import pydap-specific tools
from pydap.net import create_session
from pydap.client import get_cmr_urls

Lets get some data URLs

For this scenario, we search for remote data files (granules) from NASA's DAYMET. The following snippet uses NASA's Common Metadata Repository to search for OPeNDAP URLs

Python
daymet_ccid = "C2531982907-ORNL_CLOUD"
time_range = [dt.datetime(2014, 5, 15), dt.datetime(2024, 5, 15)] 

cmr_urls = get_cmr_urls(ccid=daymet_ccid, time_range=time_range, limit=1000)

prcp_na_urls = [url for url in cmr_urls if url.split(".nc")[0].split("Daymet_Annual_V4R1.daymet_v4_")[-1].startswith("prcp_annttl_na")]

The above code snippet should produce about 11 OPeNDAP URLs!

How to search and discover OPeNDAP URLs from NASA's Earthdata Cloud
To better understand how to search and discover OPeNDAP URLs from NASA's Earthdata Cloud, follow this link to check an entire tutorial!

Authenticate

To access data from NASA's Earthdata Cloud, you need to have an Earthdata Login (EDL) account and, importantly, a token. The next snippet of code helps with authentication, and will generate a requests session object with the necessary credentials to access data from a remote NASA server.

Python
from earthaccess.exceptions import LoginStrategyUnavailable
try:
    auth = earthaccess.login(strategy="netrc", persist=True) # you will be prompted to add your EDL credentials
except LoginStrategyUnavailable:
    auth = earthaccess.login(strategy="interactive", persist=True)

# pass Token Authorization to a new Session.
my_session = create_session(session=auth.get_session(),use_cache=True, cache_kwargs={"cache_name":'debug_case'})
my_session.cache.clear()

The snippet above shows that earthaccess was used to authenticate, which then stores the credentials in your local machine inside a stardard netrc file. Then PyDAP's create_session extracts the credentials from earthaccess to generate a requests_cache session. This is usually done to debug a session as it allows to inspects the requests done to the server! (There will be a sqlite object created in your local machine called debug.sqlite. This can be reused for later.)

At last, we create an aggregated Xarray Dataset from only 2 of the OPeNDAP URLs as follow

Python
urls = [url.replace("https","dap4") for url in prcp_na_urls]
ds = xr.open_mfdataset(urls[:2], engine="pydap", session=my_session, parallel=True, concat_dim="time", combine='nested')
ds
DAP4 or HTTPS?
The DAP4 protocol is NOT an alternative protocol for HTTPS. Replacing https with dap4 is only done to instruct Xarray (and pydap) which protocol to use. This is, how to construct the URLs to fetch data from the remote server. The URLs sent to the server will begin with https (you will verify this later).

Finally, we verify the requests sent to the server via our Cached session

Because we used PyDAP as the Xarray engine, and we used a Cached session (the my_session), we can inspect all the requests sent to the remote server

NOTE: PyDAP safely encodes the http-requests sent to the remote OPeNDAP server. Below we de-encode to show the dap4 requests made to the server

Python
[url.replace("%2F","/").replace("%5B", "[").replace("%5D","]").replace("%3A",":").replace("%3B",";") for url in my_session.cache.urls()]
Jupyter Cell Output

['https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2014.nc.dap?dap4.ce=/lambert_conformal_conic&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2014.nc.dap?dap4.ce=/lat[0:1:8074][0:1:7813]&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2014.nc.dap?dap4.ce=/lon[0:1:8074][0:1:7813]&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2014.nc.dap?dap4.ce=/time[0:1:0];/y[0:1:8074];/x[0:1:7813]&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2014.nc.dap?dap4.ce=/time_bnds[0:1:0][0:1:0]&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2014.nc.dap?dap4.ce=/time_bnds[0:1:0][1:1:1]&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2014.nc.dmr',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2015.nc.dap?dap4.ce=/lambert_conformal_conic&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2015.nc.dap?dap4.ce=/lat[0:1:8074][0:1:7813]&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2015.nc.dap?dap4.ce=/lon[0:1:8074][0:1:7813]&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2015.nc.dap?dap4.ce=/time[0:1:0];/y[0:1:8074];/x[0:1:7813]&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2015.nc.dap?dap4.ce=/time_bnds[0:1:0][0:1:0]&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2015.nc.dap?dap4.ce=/time_bnds[0:1:0][1:1:1]&dap4.checksum=true',
 'https://opendap.earthdata.nasa.gov/collections/C2531982907-ORNL_CLOUD/granules/Daymet_Annual_V4R1.daymet_v4_prcp_annttl_na_2015.nc.dmr']

NOTE: There is no dap4 scheme in the URL.

You can see from the output above that the coordinates lat and lon, each a 2D array, were downloaded from each remote file (4 extra GET requests). These are huge arrays, identical across each remote file, and downloading these coordinates can dominate the time it takes to create the dataset aggregation. There are additional extra requests (2 per remote granule) associated with the variable time_bnds.

The example shows that there are 7 requests per granule, making a total of 14 requests for only 2 remote files! This is larger than the minimum estimate of 2*(2+2)=8 GET requests (N=2, M=2, where M are the number of non-dimension coordinates, i.e. lat and lon).

This is obviously NOT a bug, but rather a safety feature built in Xarray. When attempting to aggregate datasets with coordinates like in this case, Xarray needs to check that the coordinates are identical before generating the aggregate dataset. The user may already know that the coordinates are the same, but the Xarray API makes no such assumption. As a result, it downloads the coordinates (each coordinate with its own data request) from every single remote URL.

This can clearly become a barrier for scientific discovery, since creating an aggregation for the 11 remote files alone may take up to 10x of minutes!!

What to do now? There has got to be a way!

There are useful ways to improve the performance! This is the subject for another post.