Repository: pydap/pydap (python)

Issue:
`consolidated_metadata` should better support concat_times longer than 1 element
`consolidated_metadata` pre-downloads all necessary dap responses (to be used by xarray), within the same authenticated session. It works really well, and  a `concat_dim` such as `time` can be defined. If there are N urls, then there will be N dap responses for the concat_dim (one per URL), and M dap responses, where M is the number of repeated dimensions across URL. For example, in Level 3/4 data, lon and lat are identical dimensions across all URLs, and so `M=2`.

It has become apparent that, in cases where the `time` dimension is longer that a single element, multiple `time` dap responses are downloaded per URL when using Xarray. For example, if `concat_dim`=`time`,

1. A dap response fort the whole `time` array/dimension. 
2. A dap response for first element of `time`, 
3. A dap response for the last element of `time`

This is relevant fo the MERRA-2 dataset, [see this dmr](https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dmr), time is an array of len >1. 

### Minimal example

```python
from pydap.client import consolidate_metadata

cache_session = create_session(use_cache=True)
urls = ['dap4://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4?dap4.ce=/T2M;/U2M;/V2M;/PS;/time;/lat;/lon',
        'dap4://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200124.nc4?dap4.ce=/T2M;/U2M;/V2M;/PS;/time;/lat;/lon'
         ]
# assigns cache keys, and downloads dims dap responses
# once for lat, lon, but N for time (since time differs across all N urls)
consolidate_metadata(urls, cache_session, concat_dim='time')
```

The downloaded / cache dap responses for time are:

```python
['https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dap?dap4.ce=lat%5B0%3A1%3A360%5D',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dap?dap4.ce=lon%5B0%3A1%3A575%5D',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dap?dap4.ce=time%5B0%3A1%3A23%5D',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dmr?dap4.ce=%2FT2M%3B%2FU2M%3B%2FV2M%3B%2FPS%3B%2Ftime%3B%2Flat%3B%2Flon',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200124.nc4.dap?dap4.ce=time%5B0%3A1%3A23%5D']
```

This is as expected. Now when I use xarray:
```python
ds = xr.open_mfdataset(urls, session=cache_session, engine='pydap', parallel=True, combine='nested', concat_dim='time')
```
The above command should takes ms. It takes ~10 secs. When I look at the cached urls (downloaded) I get:

```python
['https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dap?dap4.ce=lat%5B0%3A1%3A360%5D',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dap?dap4.ce=lon%5B0%3A1%3A575%5D',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dap?dap4.ce=time%5B0%3A1%3A0%5D',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dap?dap4.ce=time%5B0%3A1%3A23%5D',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dap?dap4.ce=time%5B23%3A1%3A23%5D',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200120.nc4.dmr?dap4.ce=%2FT2M%3B%2FU2M%3B%2FV2M%3B%2FPS%3B%2Ftime%3B%2Flat%3B%2Flon',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200124.nc4.dap?dap4.ce=time%5B0%3A1%3A0%5D',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200124.nc4.dap?dap4.ce=time%5B0%3A1%3A23%5D',
 'https://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200124.nc4.dap?dap4.ce=time%5B23%3A1%3A23%5D']
```

I get 2 extra dap responses per time URL. For example, the Constraint Expressions associated with these extra dap responses are:

```python
"... ?dap4.ce=time%5B0%3A1%3A0%5D",
 # and
"?dap4.ce=time%5B23%3A1%3A23%5D"
```
These look like the 1st and last element of the time dimension array. 

### What to do
This seems to be the way that `Xarray` works, to check for consistency when concatenating (in parallel) multiple datasets. For pydap, then, `consolidate_metadata` needs to check if time is longer than a single element, and if so add 2 extra dap urls for URL (first element and the last). This should trivial since there is already a length-check on the `concat_dim`.





Pull request description:
Improves support for concat_dims with size larger than 1
<!-- Thank you very much for your hard work and contribution! -->

The following Pull Request:

- [x] Closes #504 
- [x] Tests are added to demonstrate the new or fixed behavior, and all tests pass
on a local environment.

###
This PR greatly reduces the number of dap responses pydap needs to download in the case multiple (N) urls that aggregate along a dimension `concat_dim` (in most cases, `time`). 
The specific scenario is when the length of the `concat_dim` > 1. In that scenario, `3` dap responses per `concat_dim` per URL are downloaded (Xarray feature when aggregating multiple urls).

This PR makes it so that only `3` dap responses (total!) are downloaded for the `concat_dim` array (e.g. `time`). This is accomplished by forcing the download of all DMR (metadata). This is possible because of the internals of Xarray.

This code assumes that all remote data have the same dimension and identical dim shapes.

### example
Consider a year of MERRA-2 data.
```python
    len(new_urls)
>>> 366
    new_urls[:2] # check the first 2 urls
>>> ['dap4://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200101.nc4?dap4.ce=/T2M;/U2M;/V2M;/PS;/lon;/time;/lat',
 'dap4://opendap.earthdata.nasa.gov/collections/C1276812863-GES_DISC/granules/M2T1NXSLV.5.12.4%3AMERRA2_400.tavg1_2d_slv_Nx.20200102.nc4?dap4.ce=/T2M;/U2M;/V2M;/PS;/lon;/time;/lat']

## new feature for such collections.
    %% time
    consolidate_metadata(new_urls, concat_dim='time', safe_mode=True, session=cache_session)
>>> datacube has dimensions {'time[0:1:23]', 'lat[0:1:360]', 'lon[0:1:575]'} , and concat dim: `time`
CPU times: user 4.93 s, sys: 1.22 s, total: 6.15 s
Wall time: 1min 2s

### generate xarray dataset
    %%time # jupyter notebook cell magic
    ds = xr.open_mfdataset(new_urls, session=cache_session, engine='pydap', parallel=True, combine='nested', concat_dim='time')
>>> CPU times: user 2.9 s, sys: 498 ms, total: 3.39 s
Wall time: 3.02 s
```
This `O(1) secs` performance is all it will take from now on, when reusing the cache_session.
 
![Screenshot 2025-05-30 at 8 36 10 AM](https://github.com/user-attachments/assets/5a2f4229-4c74-415b-aa12-43317a3ae4c4)



## Old behavior vs New Behavior

With this new behavior the amount of URLs that are downloaded are:
```python
### new feature
    len(cache_session.cache.urls())
>>> 371
```

On the other hand, clearing the `cache_session` and skipping the `consolidate_metadata` step, the timings are:
```python
    %%time
    ds = xr.open_mfdataset(new_urls, session=cache_session, engine='pydap', parallel=True, combine='nested', concat_dim='time')
>>> CPU times: user 21 s, sys: 8.19 s, total: 29.2 s
Wall time: 30min 13s
```
and the total of downloaded URLs are:

```python
    len(cache_session.cache.urls())
>>> 2196 # !!!
```


Overall pydap can now be  10 times faster, because it downloads an order or magnitude fewer downloads





Patch:
diff --git a/src/pydap/client.py b/src/pydap/client.py
index 72775b8..3bd0a07 100644
--- a/src/pydap/client.py
+++ b/src/pydap/client.py
@@ -156,9 +156,7 @@ def open_url(
     return dataset
 
 
-def consolidate_metadata(
-    urls, session, concat_dim=None, safe_mode=False, verbose=False
-):
+def consolidate_metadata(urls, session, concat_dim=None, safe_mode=True, verbose=False):
     """Consolidates the metadata of a collection of OPeNDAP DAP4 URLs belonging to
     data cube, i.e. urls share identical variables and dimensions. This is done
     by caching the DMR response of each URL, and the DAP response of all dimensions
@@ -170,8 +168,6 @@ def consolidate_metadata(
         The URLs of the datasets that define a datacube. Each URL must begin
         with the same base URL, and begin with `dap4://`.
     session : requests-cache.CachedSession
-        A requests-cache session object. Currently, only the sqlite and memory
-        backend are fully tested. The filesystem backend is not yet supported.
     concat_dim : str, optional (default=None)
         A dimensions name (string) to concatenate across the datasets to form
         a datacube. If `None`, each dimension across the datacube
@@ -179,14 +175,16 @@ def consolidate_metadata(
         When `concat_dim` is provided, no cache key is created for that
         dimension, and the dap response associated with that dimension
         is then downloaded for each URL.
-    safe_mode : bool, optional (default=False)
-        If `True`, the DMR response is downloaded for each URL, creating a
-        empty pydap dataset. dimensions are checked for datacube consistency.
-        When `False`, only the first URL DMR response is
+    safe_mode : bool, optional (default=True)
+        If `True`, all DMR responses are downloaded for each URL, creating a
+        empty pydap dataset. dimensions names and sizes are checked for
+        datacube consistency. When `False`, only the first URL DMR response is
         downloaded, and the rest of the DMRs are assigned the same
         cache key as the first URL, to avoid downloading the DMR
-        response for each URL. This is much faster, but does not check
+        response for each URL. This is faster, but does not check
         for consistency across the URLs.
+        `NOTE`: If `concat_dim` is defined, and its dimension has a lenght
+        greater than one, `safe_mode` is automatically set to `True` always.
     verbose: bool, optional (default=False)
         For debugging purposes. If `True`, prints various URLs, normalized
         cache-keys, and other information.
@@ -223,18 +221,27 @@ def consolidate_metadata(
     dmr_urls = [
         url + ".dmr" if "?" not in url else url.replace("?", ".dmr?") for url in URLs
     ]
+    pyds = open_dmr(dmr_urls[0], session=session)
+    if concat_dim and pyds.dimensions[concat_dim] > 1:
+        if not safe_mode:
+            warnings.warn(
+                f"Length of dim `{concat_dim}` is greater than one, "
+                "reverting to `safe_mode=True`.",
+                UserWarning,
+                stacklevel=3,
+            )
+            safe_mode = True
     if safe_mode:
         with session as Session:  # Authenticate once
             with ThreadPoolExecutor(max_workers=ncores) as executor:
                 results = list(
                     executor.map(lambda url: open_dmr(url, session=Session), dmr_urls)
                 )
-        if [ds.dimensions for ds in results].count(results[0].dimensions) != len(
-            results
-        ):
+        _dim_check = results[0].dimensions
+        if not all(d == _dim_check for d in [ds.dimensions for ds in results]):
             warnings.warn(
-                "The dimensions of the datasets are not identical across all urls. "
-                "Please check the URLs and try again."
+                "The dimensions of the datasets are not identical across all remote "
+                "dataset. Please check the URLs and try again."
             )
             return None
     else:
@@ -248,21 +255,41 @@ def consolidate_metadata(
         # But needs to run so the URL is assigned the key.
         with session as Session:
             _ = download_all_urls(Session, dmr_urls, ncores=ncores)
-
     # Download dimensions once and construct cache key their dap responses
     base_url = URLs[0].split("?")[0]
     dims = set(list(results[0].dimensions.keys()))
+    add_dims = set()
     if concat_dim is not None and set([concat_dim]).issubset(dims):
         dims.remove(concat_dim)
         concat_dim_urls = [
             url.split("?")[0]
             + ".dap?dap4.ce="
             + concat_dim
-            + "[0:1:"
+            + "%5B0:1:"
             + str(results[0].dimensions[concat_dim] - 1)
-            + "]"
+            + "%5D"
             for i, url in enumerate(URLs)
         ]
+        if results[0].dimensions[concat_dim] > 1:
+            _size = results[0].dimensions[concat_dim] - 1
+            index_slices = ["%5B0:1:0%5D", f"%5B{_size}:1:{_size}%5D"]
+            concat_dim_urls += [
+                url.split("?")[0] + ".dap?dap4.ce=" + concat_dim + index
+                for url in URLs
+                for index in index_slices
+            ]
+            add_dims = set(
+                [
+                    concat_dim + "[0:1:0]",
+                    concat_dim
+                    + "["
+                    + str(results[0].dimensions[concat_dim] - 1)
+                    + ":1:"
+                    + str(results[0].dimensions[concat_dim] - 1)
+                    + "]",
+                ]
+            )
+            dims.update([concat_dim])
     else:
         concat_dim_urls = []
 
@@ -284,6 +311,7 @@ def consolidate_metadata(
     )
     if dims:
         print("datacube has dimensions", dim_ces, f", and concat dim: `{concat_dim}`")
+        dim_ces.update(add_dims)
         patch_session_for_shared_dap_cache(
             session, shared_vars=dim_ces, known_url_list=URLs, verbose=verbose
         )
@@ -582,7 +610,7 @@ def try_generate_custom_key(request, config, verbose=False):
 
     if verbose:
         print("================ request url ========================")
-        print("request.url")
+        print(request.url)
 
     shared_ext = ext
 
@@ -953,18 +981,16 @@ def get_cmr_urls(
         cmr_response["items"][i]["umm"]["RelatedUrls"]
         for i in range(len(cmr_response["items"]))
     ]
-    granule_urls = list(
-        {
-            d["URL"]
-            for item in items
-            for d in item
+    granules_urls = []
+    for item in items:
+        for i in range(len(item)):
             if (
-                d.get("Description") == "OPeNDAP request URL"
-                or d.get("Subtype") == "OPENDAP DATA"
-            )
-        }
-    )
-    return granule_urls
+                item[i].get("Description") == "OPeNDAP request URL"
+                or item[i].get("Subtype") == "OPENDAP DATA"
+            ):
+                granules_urls.append(item[i]["URL"])
+
+    return granules_urls
 
 
 if __name__ == "__main__":


Test patch:
diff --git a/src/pydap/tests/test_client.py b/src/pydap/tests/test_client.py
index e075e25..d5129d1 100644
--- a/src/pydap/tests/test_client.py
+++ b/src/pydap/tests/test_client.py
@@ -536,33 +536,51 @@ def test_cached_consolidate_metadata_inconsistent_dims(urls, safe_mode):
 )
 @pytest.mark.parametrize("concat_dim", [None, "TIME"])
 def test_consolidate_metadata_concat_dim(urls, concat_dim):
-    """Test the behavior of the chaching implemented in `consolidate_metadata`.
-    The `concat_dim` parameter means that the datasets are concatenated along
-    the specified dimension. If `concat_dim` is None, the dat response of
-    all dimensions are cached associated with that of the first URLs. This happens
-    when there is no concat dimension along all urls (sometime time is not defined),
-    and dimensions and metadata is identical across (level 3/4 datasets).
-
-    When there is a concat dimension, the dap response of all urls are all individually
-    downloaded and cached. So there should be N dap responses for N urls associated with
-    a dimension.
+    """Test the behavior of the chaching implemented in `consolidate_metadata`
+    when there is a concat dimension, and (extra) this concat_dim may be an array
+    of length >= 1.
+
+    If `concat_dim` is None, only 1 dap response per dimension is cached (1 URL),
+    and a special cache key is created for the rest of URLs. If `concat_dim` is set,
+    then N dap responses for that dimension are cached, where N is the number of URLs.
+    The rest of the (non-concat) dimensions behave the same as when `concat_dim` is
+    None.
+
+    In addition, when the `concat_dim` has length >1, there is an extra behaviour to
+    test:
+        * `safe_mode=True`: 2 extra dap responses of the `concat_dim` are downloaded,
+            per URL. These are 1st and last elements of the `concat_dim` array. Total
+            extra dap responses downloaded are 2 * N, where N is the number of URLs.
     """
     cached_session = create_session(use_cache=True, cache_kwargs={"backend": "memory"})
     cached_session.cache.clear()
-    pyds = open_dmr(urls[0].replace("dap4", "http") + ".dmr")
-    dims = list(pyds.dimensions)
-    # download all dmr for testing - not recommmended in production.
+    # download all dmr for testing - not most performant
     consolidate_metadata(
-        urls, session=cached_session, safe_mode=True, concat_dim=concat_dim
+        urls, session=cached_session, safe_mode=False, concat_dim=concat_dim
     )
+    pyds = open_dmr(urls[0].replace("dap4", "http") + ".dmr")
+    dims = list(pyds.dimensions)
+
+    N_dmr_urls = 1  # Since `safe_mode=False`, only 1 DMR is downloaded
+
     if not concat_dim:
-        assert len(cached_session.cache.urls()) == len(urls) + len(dims)
+        # Without `concat_dim` set, only one dap response is downloaded per URL.
+        assert len(cached_session.cache.urls()) == N_dmr_urls + len(dims)
     else:
-        N_concat_dims = len(urls)
+        # concat dim is set. Must download N dap responses for the concat_dim.
+        N_concat_dims = len(urls)  # see below !
+        if pyds.dimensions[concat_dim] > 1:
+            N_concat_dims = 1  # only one dap response of the concat dim is downloaded!
+
+            # there are 2 extra dap responses: 1st and last elements
+            # of concat dim array (one url only!)
+            N_concat_dims += 2
+            N_dmr_urls = len(urls)  # safe_mode becomes True, and all dmrs
+            # are downloaded
         N_non_concat_dims = len(dims) - 1  # only one dimension is concatenated
         assert (
             len(cached_session.cache.urls())
-            == len(urls) + N_concat_dims + N_non_concat_dims
+            == N_dmr_urls + N_concat_dims + N_non_concat_dims
         )