Hyrax 1.18.0 is now available.
Read the release notes →
Scalable OPeNDAP access with Xarray: Avoiding full-array downloads when subsetting
Xarray is an incredibly useful library that hides the complexity of accessing, subsetting, plotting, and streaming data behind its simple and intuitive methods. It can also read data through pluggable backend engines; specifying for example engine="pydap", lets users access remote data as if they were local files. However, this convenience comes at a cost: adding multiple libraries into a single toolchain introduces inefficiencies that are not always visible to the user and hard to diagnose. Here, I present a workflow that avoids a hidden performance bottleneck when accessing aggregated data with Xarray and OPeNDAP. Because Xarray’s default behavior requests the entire array before subsetting, anyone who works with remote data as if it were local will encounter this bottleneck. The fix requires only one additional line, which ensures that subsetting is performed remotely by the OPeNDAP server.
By default, Xarray requests the entire array in an aggregation, and then subsets locally
Xarray treats each remote OPeNDAP URL as a single chunk. Because Xarray works chunk by chunk, it downloads the full array from each file in the aggregation and only then subsets it locally, instead of requesting just the subset from the server.
To avoid this behavior, we need to specify chunks that align with the subset when opening the dataset.
When accessing remote, distributed data, performance depends largely on how much data is transferred (see Figure 1a). Ideally, a client requests only the values it needs. With Xarray’s default behavior, however, each request can be far larger than necessary, limiting scalability.
Optimized streaming workflow
Below, we present a simple modification to a common Xarray workflow that ensures subsetting happens on the remote OPeNDAP server before the data are sent. It adds one additional step: prechunk when opening the Xarray Dataset. The resulting workflow is represented in Figure 1. This small change substantially improves the scalability of data access.
Python
#===================================================
# Identify the size of subset (from coordinate data)
#===================================================
# download entire coordinate first
Lat, Lon = ...
# size of entire coordinate arrays
nlat = len(Lat)
nlon = len(Lon)
# Create masking arrays
iLat = np.where(...) # 1D array
iLon = np.where(...) # 1D array
# create chunk mapping, making sure that it coincides with the chunk boundaries
Chunks={
"lat": (iLat[0], len(iLat), nlat - iLat[0] - len(iLat)),
"lon": (iLon[0], len(iLon), nlon - iLon[0] - len(iLon)),
}
# (Pre-chunk) when creating aggregate Xarray dataset
ds = xr.open_mfdataset(
urls,
engine="pydap",
chunks=Chunks # < ------------- see here
parallel=True,
...,
)
Specifying chunks when opening the dataset (e.g., xr.open_dataset(url, engine="pydap", chunks=...)), with chunk sizes and boundaries that match the expected subset (Fig. 1a), sets Xarray to issue a single data request per file.
NOTE
Pre-chunking differs from re-chunking an already opened dataset with ds.chunk(), which does not change how the data are requested from the server.
Below we demonstrate with data and a reproducible example
1. Create an isolated environment and launch a jupyter notebook
Terminal
$ mamba create -n opendap_env -c conda-forge python=3.12 ipython pydap jupyterlab earthaccess netCDF4 xarray dask
$ mamba activate opendap_env
$ jupyter lab2. Import necessary libraries
Python
import xarray as xr
import datetime as dt
import numpy as np
import earthaccess
# import pydap-specific tools
from pydap.net import create_session, open_url
from pydap.client import get_cmr_urls
3. Discover OPeNDAP URLs to access temperature reanalysis data from MERRA2
Python
MERRA2_M2I3NPASM_ccid = "C1276812879-GES_DISC" # short name M2I3NPASM
# Fall data from 2 years
time_ranges = [[dt.datetime(year, 9, 21), dt.datetime(year, 12, 22)] for year in range(2024, 2026)]
cmr_urls = [urls for time in time_ranges for urls in get_cmr_urls(ccid=MERRA2_M2I3NPASM_ccid, time_range=time, limit=1000)]
There are 186 URLs found. For this demonstration we will only use 10.
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!
4. Authenticate via EDL and earthaccess
Before downloading data, we need to authenticate. For that, you must have a valid EDL account and have created a token.
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()
Having authenticated and passed along the credentials to a requests-cache session object, we can now begin to use Xarray to request data from the OPeNDAP server.
5. Use PyDAP to inspect the remote datafile and accelerate the aggregate dataset creation
The steps below are a set of additional optimizations regarding the creation of an aggregate Xarray Dataset, when the data files are remote. Tthe optimization steps are discussed in the previous post, but we describe them in the code block below.
Python
# Inspect a single remote file
url = cmr_urls1[0].replace("https","dap4")
pyds = open_url(url, session=my_session)
# ========
# 1) Identify variable of interest, and its dimensions
Variables = ["/T"]
dims = list(set([dim for var in Variables for dim in pyds[var[1:]].dims]))
# ========
# 2) Keep all dimensions of Temperature in first URL
url1 = cmr_urls[0].replace("https","dap4") + "?dap4.ce=" + ";".join(dims+Variables)
# ========
# 3) For rest of URLs, drop spatial dimensions.
# spatial dims are all identical across files.
urls_rest = [url.replace("https","dap4")+"?dap4.ce=" + ";".join(["/time"]+Variables) for url in cmr_urls[1:10]]
# ======
# 4) Final list of OPeNDAP URLs
urls = [url1] + urls_rest
6. Identify a subset.
This requires creating an Xarray dataset. We do not need to aggregate yet. Because this is the case of Level 3 data, the aggregation will form a datacube, and any slicing will be applied equally across all files.
Consider the scenario where we only want a small region in the North Atlantic Ocean, bounded by the spatial range: -75 < longitude < 15, and 15 < latitude < 65. The following snippet of code identifies the slice that streams the desired subset to a new file on disk.
Python
# Look at the 1st URL which has all dimension data
ds1 = xr.open_dataset(url1, engine="pydap")
#### =====================================
# Identify subset given the lat/lon bounds
#### =====================================
# Min/max of lon values & Min/Max of lat values
lat_min, lat_max = 15, 65
lon_min, lon_max = -75, 15
Lon, Lat = ds1['lon'], ds1['lat']
iLon = np.where((Lon>lon_min)&(Lon < lon_max))[0]
iLat = np.where((Lat>lat_min)&(Lat < lat_max))[0]
7. Create the aggregated Xarray dataset: Recommended workflow
Python
# clear any cache data
my_session.cache.clear()
# Need to chunk according to the size and location of subset
nlat = len(Lat)
nlon = len(Lon)
Chunks={
"lat": (iLat[0], len(iLat), nlat - iLat[0] - len(iLat)),
"lon": (iLon[0], len(iLon), nlon - iLon[0] - len(iLon)),
}
# Create the aggregate Xarray Dataset with an extra chunking parameter to match the subset
ds = xr.open_mfdataset(
urls, engine="pydap", session=my_session,
concat_dim='time', parallel=True, combine='nested',
chunks=Chunks, # <--------- ! Must add this
)
# now try again streaming the data from OPeNDAP servers
ds.isel(lon=slice(iLon[0], iLon[-1]+1), lat=slice(iLat[0], iLat[-1]+1)).to_netcdf("data/test2.nc4")
We now inspect the data responses that were requested from the remote OPeNDAP server.
Python
[url.replace("%2F", "/").replace("%5B","[").replace("%5D","]").replace("%3A",":") for url in my_session.cache.urls()]
Jupyter Cell Output
['https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240921.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][211:1:309][169:1:311]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240922.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][211:1:309][169:1:311]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240923.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][211:1:309][169:1:311]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240924.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][211:1:309][169:1:311]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240925.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][211:1:309][169:1:311]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240926.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][211:1:309][169:1:311]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240927.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][211:1:309][169:1:311]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240928.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][211:1:309][169:1:311]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240929.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][211:1:309][169:1:311]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240930.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][211:1:309][169:1:311]&dap4.checksum=true']
The query part of the URL (anything after a ?) specifies the data request.
Where is the subset taking place?
Our approach ensures a single data request is made, so that subsetting always takes place close to the data!
What if you do not pre-chunk?
This is the default behavior. We show below that the entire array is requested
Python
# clear the cached data
my_session.cache.clear()
# create dataset aggregate (no prechunking)
ds = xr.open_mfdataset(
urls, engine="pydap",
session=my_session, concat_dim='time',
parallel=True, combine='nested',
)
ds.isel(lon=slice(iLon[0], iLon[-1]+1), lat=slice(iLat[0], iLat[-1]+1)).to_netcdf("data/test2.nc4")
Inspecting the URL as we did before, we get
Jupyter Cell Output
['https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240921.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][0:1:360][0:1:575]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240922.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][0:1:360][0:1:575]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240923.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][0:1:360][0:1:575]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240924.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][0:1:360][0:1:575]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240925.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][0:1:360][0:1:575]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240926.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][0:1:360][0:1:575]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240927.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][0:1:360][0:1:575]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240928.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][0:1:360][0:1:575]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240929.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][0:1:360][0:1:575]&dap4.checksum=true',
'https://opendap.earthdata.nasa.gov/collections/C1276812879-GES_DISC/granules/M2I3NPASM.5.12.4:MERRA2_400.inst3_3d_asm_Np.20240930.nc4.dap?dap4.ce=/T[0:1:7][0:1:41][0:1:360][0:1:575]&dap4.checksum=true']
We can see that the entire array was downloaded. This means that all the subsetting took place locally. This represents a bottleneck to scalability.
