Repository: pydap/pydap
Pull request #428 (2024-12-20 00:57:38)

Decode chunk header to set length of DMR in a generic way and endianness of data
<!-- Thank you very much for your hard work and contribution! -->

The following Pull Request: ~(DRAFT)~

- [x] Closes an already existing opened issue #xxxx for fixing a bug, requesting
a new feature, etc.
- [x] Tests are added to demonstrate the new or fixed behavior, and all tests pass
on a local environment.

Splits the dap response into dmr (first chunk) and data in a generic way, taking into account the first 4 bytes of the response. 


Mikejmnez: with these changes, pydap can set endianness correctly only once (from 1st chunk header), and can correctly read basic DAP4 data from tds and hyrax. Reading from these two

```python
url_tds = "https://thredds-test.unidata.ucar.edu/thredds/dap4/dev/d4icomp/SimpleGroup.nc4.dap?dap4.ce=/SimpleGroup/X&dap4.checksum=true"
url_hds = 'http://test.opendap.org:8080/opendap/dap4/SimpleGroup.nc4.h5.dap?dap4.ce=/SimpleGroup/X'
```

X is an array on int16 spanning from 0 to 39.

With pydap I get:
```python
dataset_tds['/SimpleGroup/X']
>>> <BaseType with data array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37, 38, 39], dtype=int16)>
```
and
```python
dataset_hds['/SimpleGroup/X']
>>> <BaseType with data array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37, 38, 39], dtype=int16)>
```


Still draft mode. I need to refactor several things. And test against other data types in that same file.
Mikejmnez: ready for review! It sets the endianness of the data by reading the first chunk (before the DMR). It makes the assumption that the endianness does not change (tds approach). I check the datasets on both tds and hyrax servers and they look identical to the original dataset (which I created). Both 1D and 2D arrays (int16 and float32).



@ndp-opendap 

Patch:
diff --git a/src/pydap/client.py b/src/pydap/client.py
index 6ebf5f5..fc3f435 100644
--- a/src/pydap/client.py
+++ b/src/pydap/client.py
@@ -46,7 +46,6 @@ lazy mechanism for function call, supporting any function. Eg, to call the
 
 from io import BytesIO, open
 
-import numpy
 from requests.utils import urlparse, urlunparse
 
 import pydap.handlers.dap
@@ -56,6 +55,7 @@ import pydap.net
 import pydap.parsers.das
 import pydap.parsers.dds
 import pydap.parsers.dmr
+from pydap.handlers.dap import UNPACKDAP4DATA
 from pydap.lib import DEFAULT_TIMEOUT as DEFAULT_TIMEOUT
 
 
@@ -105,51 +105,24 @@ def open_file(file_path, das_path=None):
         return open_dmr_file(file_path=file_path)
 
 
-def get_dmr_length(file_path):
-    with open(file_path, "rb") as f:
-        # First two bytes are CRLF
-        if f.peek()[0:2] == b"\x04\x00":
-            f.seek(2)
-            dmr_len = numpy.frombuffer(f.read(2), dtype=">u2")[0]
-        else:
-            dap = f.read()
-            dmr = b"<?xml" + dap.split(b"<?xml")[1]
-            dmr = dmr.split(b"</Dataset>")[0] + b"</Dataset>\n\r\n"
-            dmr_len = len(dmr)
-    return dmr_len
-
-
 def open_dmr_file(file_path):
-    dmr_len = get_dmr_length(file_path)
+    """
+    Opens a DMR. This differs from a dap response, since it is a single chunk,
+    and there is no chunk header at the top of the file.
+    """
     with open(file_path, "rb") as f:
-        if f.peek()[0:2] == b"\x04\x00":
-            # First 2 bytes are CRLF, second two bytes give the length of the
-            # DMR; we skip over them
-            f.seek(4)
-            # We read the DMR minus the CRLF and newline (3 bytes)
-        dmr = f.read(dmr_len)
+        dmr = f.read()
     dmr = dmr.decode("ascii")
     dataset = pydap.parsers.dmr.dmr_to_dataset(dmr)
     return dataset
 
 
-def open_dap_file(file_path, das_path=None):
-    """Open a file downloaded from a `.dap` (dap4) response, retunring a
-    dataset.
-
-    Optionally, read also the `.das` response to assign attributes to the
+def open_dap_file(file_path):
+    """Open a file downloaded from a `.dap` (dap4) response, returning a
     dataset.
     """
-
-    dataset = open_dmr_file(file_path)
-
     with open(file_path, "rb") as f:
-        dmr_len = get_dmr_length(file_path)
-        # if f.peek()[0:2] == b'\x04\x00':
-        f.seek(dmr_len)
-        _ = f.read(4)
-        pydap.handlers.dap.unpack_dap4_data(f.read(), dataset)
-    return dataset
+        return UNPACKDAP4DATA(f).dataset
 
 
 def open_dods_file(file_path, das_path=None):
diff --git a/src/pydap/handlers/dap.py b/src/pydap/handlers/dap.py
index 4475dbe..8855739 100644
--- a/src/pydap/handlers/dap.py
+++ b/src/pydap/handlers/dap.py
@@ -19,11 +19,12 @@ import sys
 import warnings
 import warnings as _warnings
 from concurrent.futures import ThreadPoolExecutor
-from io import BytesIO
+from io import BufferedReader, BytesIO
 from itertools import chain
 
 import numpy
 from requests.utils import urlparse, urlunparse
+from webob.response import Response
 
 from pydap.handlers.lib import BaseHandler, ConstraintExpression, IterData
 from pydap.lib import (
@@ -313,42 +314,6 @@ def safe_dds_and_data(r, user_charset):
     return dds.decode(get_charset(r, user_charset)), data
 
 
-def safe_dmr_and_data(r, user_charset, url):
-    if r.content_encoding == "gzip":
-        raw = gzip.GzipFile(fileobj=BytesIO(r.body)).read()
-    else:
-        raw = r.body
-    logger.info("Saving and splitting dmr+")
-    try:
-        dmr, data = raw.split(b"</Dataset>", 1)
-    except ValueError:
-        logger.exception("Failed to split the following DMR+ \n %s" % raw)
-        import codecs
-        import pickle
-
-        picked_response = str(codecs.encode(pickle.dumps(r), "base64").decode())
-        header = "pickled response (base64)"
-        logger.exception(
-            header
-            + ": \n ----BEGIN PICKLE----- \n %s \n -----END PICKLE-----"
-            % picked_response
-        )
-
-    dmr = dmr[4:] + b"</Dataset>"
-    dmr = dmr.decode(get_charset(r, user_charset))
-    if "thredds" in url.split("/") or "dap4" in url.split("/"):
-        data = data[2:]  # TDS
-        warnings.warn(
-            "Full DAP4 support for TDS data servers remains under development"
-            " for pydap. We recommend to continue to use DAP2 protocol instead."
-            " See: https://github.com/pydap/pydap/issues/401",
-            UserWarning,
-        )
-    else:
-        data = data[3:]  # Hyrax
-    return dmr, data
-
-
 class BaseProxyDap2(object):
     """A proxy for remote base types.
 
@@ -495,12 +460,7 @@ class BaseProxyDap4(BaseProxyDap2):
         )
 
         raise_for_status(r)
-        dmr, data = safe_dmr_and_data(r, self.user_charset, self.baseurl)
-
-        # Parse received dataset:
-        dataset = dmr_to_dataset(dmr)
-        dataset = unpack_dap4_data(data, dataset)
-
+        dataset = UNPACKDAP4DATA(r, self.user_charset).dataset
         self.checksum = dataset[self.id].attributes["checksum"]
         self.data = dataset[self.id].data
         return self.data
@@ -784,7 +744,39 @@ def unpack_dap2_data(xdr_stream, dataset):
     return unpack_children(xdr_stream, dataset)
 
 
+def find_pattern_in_string_iter(pattern, i):
+    last_chunk = b""
+    length = len(pattern)
+    for this_chunk in i:
+        last_chunk += this_chunk
+        m = re.search(pattern, last_chunk)
+        if m:
+            return last_chunk[m.end() :]
+        last_chunk = last_chunk[-length:]
+
+
+def dump():  # pragma: no cover
+    """Unpack dods response into lists.
+
+    Return pretty-printed data.
+
+    """
+    dods = sys.stdin.read()
+    dds, xdrdata = dods.split(b"\nData:\n", 1)
+    dataset = dds_to_dataset(dds)
+    xdr_stream = io.BytesIO(xdrdata)
+    data = unpack_dap2_data(xdr_stream, dataset)
+    pprint.pprint(data)
+
+
 def decode_chunktype(chunk_type):
+    """
+    Takes the chunk_type from chunk header embeded in the dap response,
+    and returns a tuple of
+        last_chunk: bool
+        error: bool
+        endian
+    """
     encoding = "{0:03b}".format(chunk_type)
     if sys.byteorder == "little":
         # If our machine's byteorder is little,
@@ -856,60 +848,98 @@ def stream2bytearray(data):
     return buffer
 
 
-def get_endianness(xdr_stream):
-    chunk_header = xdr_stream.peek(4)[0:4]
+def get_endianness(chunk_header):
     chunk_header = numpy.frombuffer(chunk_header, dtype=">u4")[0]
     chunk_type = (chunk_header >> 24) & 0xFF
-    last, error, endian = decode_chunktype(chunk_type)
+    _, _, endian = decode_chunktype(chunk_type)
     return endian
 
 
-def unpack_dap4_data(data, dataset):
-    endian = get_endianness(BytesReader(data))
-    checksum_dtype = numpy.dtype(endian + "u4")
-    buffer = stream2bytearray(data)
-    start = 0
-    for variable in walk(dataset, BaseType):
-        # variable_name = variable.name
-        # variable = dataset[variable_name]
-        count = get_count(variable)
-        stop = start + count
-        data = decode_variable(
-            buffer, start=start, stop=stop, variable=variable, endian=endian
-        )
-        checksum = numpy.frombuffer(
-            buffer[stop : stop + 4], dtype=checksum_dtype
-        ).byteswap("=")
-        if isinstance(variable, BaseType):
-            variable._set_data(data)
-        elif isinstance(variable, GridType):
-            variable._set_data([data.data])
-        variable.attributes["checksum"] = checksum
-        # Jump over the 4 byte chunk_header
-        start = stop + 4
-    return dataset
-
-
-def find_pattern_in_string_iter(pattern, i):
-    last_chunk = b""
-    length = len(pattern)
-    for this_chunk in i:
-        last_chunk += this_chunk
-        m = re.search(pattern, last_chunk)
-        if m:
-            return last_chunk[m.end() :]
-        last_chunk = last_chunk[-length:]
-
-
-def dump():  # pragma: no cover
-    """Unpack dods response into lists.
-
-    Return pretty-printed data.
-
+class UNPACKDAP4DATA(object):
     """
-    dods = sys.stdin.read()
-    dds, xdrdata = dods.split(b"\nData:\n", 1)
-    dataset = dds_to_dataset(dds)
-    xdr_stream = io.BytesIO(xdrdata)
-    data = unpack_dap2_data(xdr_stream, dataset)
-    pprint.pprint(data)
+    Unpacks DAP4 response, remote or local, which is split into chunks. The
+    first chunk contains the DMR response, and the endianness is defined in the first
+    4 bytes before the DMR. Once the endianness is set, it remains unchanged. This
+    makes the assumption that the all variables within dataset have the same endianness.
+
+    Parameters:
+    -----------
+        r: contains the dap response.
+            May be a Webob.response.Response created from pydap.net.GET if the dataset
+            is remote (from a url), or a `io.BufferedReader` if the data is local within
+            a filesystem. See `pydap.net.get.open_dap_file`
+    """
+
+    def __init__(self, r, user_charset="ascii"):
+        self.user_charset = user_charset
+        if isinstance(r, Response):  # a Webob response
+            self.r = r
+            if self.r.content_encoding == "gzip":
+                self.raw = BytesReader(
+                    gzip.GzipFile(fileobj=BytesIO(self.r.body)).read()
+                )
+            else:
+                self.raw = BytesReader(r.body)
+        elif isinstance(r, BufferedReader):
+            # r comes from reading a local file
+            self.r = Response()  # make empty response
+            self.raw = BytesReader(r.read())
+        else:
+            raise TypeError(
+                """
+                Unrecognized file type object for unpacking dap4 binary data.
+                Acceptable formats are `webob.response.Response` and
+                `io.BufferedReader`
+                """
+            )
+        self.dmr, self.data, self.endianness = self.safe_dmr_and_data()
+        # need to split dmr from data
+        dataset = dmr_to_dataset(self.dmr)
+        self.dataset = self.unpack_dap4_data(dataset)
+
+    def safe_dmr_and_data(self):
+        """
+        Splits the dap response (.dap) into the dmr (metadata), and the raw
+        (binary) data. It also computes the endianness of the data.
+
+        Returns:
+            dmr, data, endianness
+        """
+        # decode the first 4 bytes are CRLF
+        chunk_header = numpy.frombuffer(self.raw.read(4), dtype=">u4")[0]
+        dmr_length = chunk_header & 0x00FFFFFF
+        chunk_type = (chunk_header >> 24) & 0xFF
+        dmr = self.raw.read(dmr_length).decode(get_charset(self.r, self.user_charset))
+        data = self.raw.data
+        # get endianness from first chunk
+        _, _, endianness = decode_chunktype(chunk_type)
+        return dmr, data, endianness
+
+    def unpack_dap4_data(self, dataset):
+        """
+        Takes a pydap.DatasetType previously created, and populates its variables
+        (BaseType only) with data that is currently in binary form (within a dap
+        response).
+        """
+        # need self. data and self.dataset
+        checksum_dtype = numpy.dtype(self.endianness + "u4")
+        buffer = stream2bytearray(self.data)
+        start = 0
+        for variable in walk(dataset, BaseType):
+            count = get_count(variable)
+            stop = start + count
+            data = decode_variable(
+                buffer,
+                start=start,
+                stop=stop,
+                variable=variable,
+                endian=self.endianness,
+            )
+            checksum = numpy.frombuffer(
+                buffer[stop : stop + 4], dtype=checksum_dtype
+            ).byteswap("=")
+            variable._set_data(data)
+            variable.attributes["checksum"] = checksum
+            # Jump over the 4 byte chunk_header
+            start = stop + 4
+        return dataset


Test patch:
diff --git a/src/pydap/tests/test_handlers_dap.py b/src/pydap/tests/test_handlers_dap.py
index ed2372a..528586f 100644
--- a/src/pydap/tests/test_handlers_dap.py
+++ b/src/pydap/tests/test_handlers_dap.py
@@ -1,14 +1,18 @@
 """Test the DAP handler, which forms the core of the client."""
 
+import os
 import unittest
 
 import numpy as np
+from webob.response import Response
 
 import pydap.model
 from pydap.handlers.dap import (
+    UNPACKDAP4DATA,
     BaseProxyDap2,
     DAPHandler,
     SequenceProxy,
+    dmr_to_dataset,
     find_pattern_in_string_iter,
     walk,
 )
@@ -600,3 +604,27 @@ class TestArrayStringBaseType(unittest.TestCase):
     def test_getitem(self):
         """Test the ``__getitem__`` method."""
         np.testing.assert_array_equal(self.data[:], self.original_data)
+
+
+class TestUnpackDap4Data(unittest.TestCase):
+    """
+    Tests some properties of the class UNPACKDAP4DATA.
+    """
+
+    def setUp(self):
+        file_path = "data/daps/coads_climatology.nc.dap"
+        abs_path = os.path.join(os.path.dirname(__file__), file_path)
+        with open(abs_path, "rb") as f:
+            self.unpacker = UNPACKDAP4DATA(f)
+
+    def testEndian(self):
+        """test the endianness of the first chunk encoded within the dap file"""
+        endian = self.unpacker.endianness
+        self.assertEqual(endian, "<")
+
+    def testDMR(self):
+        ds = dmr_to_dataset(self.unpacker.dmr)
+        self.assertEqual(ds.variables(), {"SST": np.dtype(">f4")})
+
+    def testResponse(self):
+        self.assertIsInstance(self.unpacker.r, Response)
diff --git a/src/pydap/tests/test_parsers_dap.py b/src/pydap/tests/test_parsers_dap.py
index a8b406c..0224a41 100644
--- a/src/pydap/tests/test_parsers_dap.py
+++ b/src/pydap/tests/test_parsers_dap.py
@@ -2,13 +2,12 @@ import os
 
 import numpy
 
-import pydap.client
+from ..client import open_dap_file
 
 
 def load_dap(file_path):
     abs_path = os.path.join(os.path.dirname(__file__), file_path)
-    dataset = pydap.client.open_dap_file(abs_path)
-    return dataset
+    return open_dap_file(abs_path)
 
 
 def test_jpl1():
diff --git a/src/pydap/tests/test_parsers_dmr.py b/src/pydap/tests/test_parsers_dmr.py
index 070c18c..9f8840d 100644
--- a/src/pydap/tests/test_parsers_dmr.py
+++ b/src/pydap/tests/test_parsers_dmr.py
@@ -7,6 +7,7 @@ from xml.etree import ElementTree as ET
 
 import numpy as np
 
+from ..client import open_dmr_file
 from ..lib import walk
 from ..model import BaseType
 from ..parsers.dmr import dmr_to_dataset, get_groups
@@ -14,10 +15,7 @@ from ..parsers.dmr import dmr_to_dataset, get_groups
 
 def load_dmr_file(file_path):
     abs_path = os.path.join(os.path.dirname(__file__), file_path)
-    with open(abs_path, "r") as dmr_file:
-        dmr = dmr_file.read()
-    dataset = dmr_to_dataset(dmr)
-    return dataset
+    return open_dmr_file(abs_path)
 
 
 class TestDMRParser(unittest.TestCase):