Repository: pydap/pydap

Issue:
`pydap` appears to be incompatible with `numpy>=2`
<!-- Thank you very much for your input! -->

### Description

Importing `pydap.client` fails if `numpy>=2` is installed (noticed in `xarray`'s CI).

### Minimal example that reproduces the bug, or a new desired feature

With `numpy>=2` in the environment:

```pytb
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 1
----> 1 import pydap.client

File .../lib/python3.11/site-packages/pydap/client.py:53
     51 import pydap.lib
     52 import pydap.net
---> 53 import pydap.handlers.dap
     54 import pydap.parsers.dds
     55 import pydap.parsers.dmr

File .../lib/python3.11/site-packages/pydap/handlers/dap.py:16
     14 import re
     15 from itertools import chain
---> 16 from numpy.lib.arrayterator import Arrayterator
     18 # handlers should be set by the application
     19 # http://docs.python.org/2/howto/logging.html#configuring-logging-for-a-library
     20 import logging

ModuleNotFoundError: No module named 'numpy.lib.arrayterator'
```

<details><summary>Environment</summary>

```
python : 3.11.9 (conda-forge)
numpy: 2.0.0
pydap: 3.4.0
```

</details>


Discussion:
Mikejmnez: With Numpy = 2, a sequence now displays data slightly different. Consider the following example:
## Setup of example
```python
        >>> import numpy as np
        >>> data = np.array([
        ... (10, 15.2, 'Diamond_St'),
        ... (11, 13.1, 'Blacktail_Loop'),
        ... (12, 13.3, 'Platinum_St'),
        ... (13, 12.1, 'Kodiak_Trail')],
        ... dtype=np.dtype([
        ... ('index', np.int32), ('temperature', np.float32),
        ... ('site', np.dtype('|S14'))]))
        ...
        >>> seq = SequenceType('example')
        >>> seq['index'] = BaseType('index')
        >>> seq['temperature'] = BaseType('temperature')
        >>> seq['site'] = BaseType('site')
        >>> seq.data = data
```
###    Iteraring over the sequence to print data:

```python
        >>> for line in seq.iterdata():
        ...     print(line)
        (np.int32(10), np.float32(15.2), 'Diamond_St')
        (np.int32(11), np.float32(13.1), 'Blacktail_Loop')
        (np.int32(12), np.float32(13.3), 'Platinum_St')
        (np.int32(13), np.float32(12.1), 'Kodiak_Trail')

```
Note the `np.type` for each element of the array. 

### Before , the same code would print

```python
        (10, 15, 'Diamond_St')
        (11, 13.1, 'Blacktail_Loop')
        (12, 13.3, 'Platinum_St')
        (13, 12.1, 'Kodiak_Trail')
```
## Why?
With numpy 2 the following code throws an error

The change alone  introduced many errors. which required the change here
https://github.com/pydap/pydap/blob/e133372aea4f033b950c1137953dfd8b97c81b88/src/pydap/responses/dods.py#L45-L46

The error arises because newbytesorder is no longer a method of `x.astype` ...



Mikejmnez: I found this simple way to decode the representation of numpy types into native Python types, with a slight modification to the following function. This would allow to revert all the changes to the documentation. That is, change the following function:
https://github.com/pydap/pydap/blob/e133372aea4f033b950c1137953dfd8b97c81b88/src/pydap/lib.py#L330-L335

To

```python
def decode_np_strings(numpy_var): 
    """Given a fixed-width numpy string, decode it to a unicode type""" 
    if isinstance(numpy_var, bytes) and hasattr(numpy_var, "tobytes"): 
        return numpy_var.tobytes().decode("utf-8") 
    else:
        if isinstance(numpy_var, np.generic):
            # starting with numpy 2.0 need to explicitly convert
            # numpy scalar to native Python type
            return numpy_var.item()
        else: 
            return numpy_var
```

However the new decode_np_strings yields now

```python
        >>> for line in seq.iterdata():
        ...     print(line)
    (10, 15.199999809265137, 'Diamond_St')
    (11, 13.100000381469727, 'Blacktail_Loop')
    (12, 13.300000190734863, 'Platinum_St')
    (13, 12.100000381469727, 'Kodiak_Trail')
```
instead of 
>
>        (10, 15, 'Diamond_St')
>        (11, 13.1, 'Blacktail_Loop')
>        (12, 13.3, 'Platinum_St')
>        (13, 12.1, 'Kodiak_Trail')
>
>

Nothing wrong - since this is the inherent precision of `np.float32`, but it may lead to the wrong interpretation that pydap is somehow changing the data.




### Thus
I decided NOT to change `np.decode_string` and instead change the documentation

Mikejmnez: @jgallagher59701 ready for merger 
Mikejmnez: Awesome, thanks!

Patch:
diff --git a/ci/environment.yml b/ci/environment.yml
index ae7046f..b85544a 100644
--- a/ci/environment.yml
+++ b/ci/environment.yml
@@ -2,7 +2,7 @@ name: pydap_tests
 channels:
 - conda-forge
 dependencies:
-- numpy < 2.0
+- numpy
 - netCDF4
 - pytest
 - pytest-cov
diff --git a/docs/developer_data_model.rst b/docs/developer_data_model.rst
index cf8d56e..35dee1b 100644
--- a/docs/developer_data_model.rst
+++ b/docs/developer_data_model.rst
@@ -288,9 +288,9 @@ We can also iterate over the ``SequenceType``. In this case, it will return a se
 
     >>> for record in s.iterdata():
     ...     print(record)
-    (1, 10)
-    (2, 20)
-    (3, 30)
+    (np.int32(1), np.int16(10))
+    (np.int32(2), np.int16(20))
+    (np.int32(3), np.int16(30))
 
 Prior to pydap 3.2.2, this approach was not possible and one had to iterate directly over ``SequenceType``:
 
@@ -298,9 +298,9 @@ Prior to pydap 3.2.2, this approach was not possible and one had to iterate dire
 
     >>> for record in s:
     ...     print(record)
-    (1, 10)
-    (2, 20)
-    (3, 30)
+    (np.int32(1), np.int16(10))
+    (np.int32(2), np.int16(20))
+    (np.int32(3), np.int16(30))
 
 This approach will be deprecated in pydap 3.4.
 
@@ -310,7 +310,7 @@ Numpy, since we can reference them by column (``s['a']``) or by index:
 .. doctest::
 
     >>> s[1].data
-    (2, 20)
+    np.void((2, 20), dtype=[('a', '<i4'), ('long%20%26%20complicated', '<i2')])
     >>> s[ s.a < 3 ].data
     array([(1, 10), (2, 20)],
           dtype=[('a', '<i4'), ('long%20%26%20complicated', '<i2')])
@@ -353,7 +353,7 @@ And here is an example of how we would use it:
            [10, 20]])>
     >>> for record in s2.iterdata():
     ...     print(record)
-    (10, 20)
+    (np.int64(10), np.int64(20))
 
 One can also iterate directly over the ``IterData`` object to obtain the data:
 
@@ -361,9 +361,11 @@ One can also iterate directly over the ``IterData`` object to obtain the data:
 
     >>> for record in s2:
     ...     print(record)
-    (10, 20)
+    (np.int64(10), np.int64(20))
 
-This approach will not be deprecated in pydap 3.4.
+This approach will not be deprecated in pydap 3.4. NOTE: For numpy > 2.0, iterating over the
+``IterData`` object returns a record specifying the individual types of the elements of the
+sequence
 
 There are many implementations of classes derived from ``IterData``: ``pydap.handlers.dap.SequenceProxy`` is a proxy to
 sequential data on Opendap servers, ``pydap.handlers.csv.CSVProxy`` wraps a CSV file,
diff --git a/src/pydap/handlers/lib.py b/src/pydap/handlers/lib.py
index ca7f542..c56827b 100644
--- a/src/pydap/handlers/lib.py
+++ b/src/pydap/handlers/lib.py
@@ -17,7 +17,7 @@ import sys
 
 import numpy as np
 import pkg_resources
-from numpy.lib.arrayterator import Arrayterator
+from numpy.lib import Arrayterator
 from webob import Request
 
 from pydap.exceptions import ConstraintExpressionError, ExtensionNotSupportedError
diff --git a/src/pydap/model.py b/src/pydap/model.py
index f706cb4..e252df2 100644
--- a/src/pydap/model.py
+++ b/src/pydap/model.py
@@ -128,15 +128,15 @@ obects are very flexible. Data can be accessed by iterating over the object::
 
     >>> for record in cast.iterdata():
     ...     print(record)
-    (10.0, 17.0, 35.0, '1')
-    (20.0, 15.0, 35.0, '2')
+    (np.float32(10.0), np.float32(17.0), np.float32(35.0), '1')
+    (np.float32(20.0), np.float32(15.0), np.float32(35.0), '2')
 
 It is possible to select only a few variables::
 
     >>> for record in cast['salinity', 'depth'].iterdata():
     ...     print(record)
-    (35.0, 10.0)
-    (35.0, 20.0)
+    (np.float32(35.0), np.float32(10.0))
+    (np.float32(35.0), np.float32(20.0))
 
     >>> cast['temperature'].dtype
     dtype('float32')
@@ -155,13 +155,13 @@ When constrained, it yields the SequenceType:
     <class 'pydap.model.SequenceType'>
     >>> for record in cast[ cast['temperature'] < 16 ].iterdata():
     ...     print(record)
-    (20.0, 15.0, 35.0, '2')
+    (np.float32(20.0), np.float32(15.0), np.float32(35.0), '2')
 
 As mentioned earlier, it is still possible to iterate directly over data::
 
     >>> for record in cast[ cast['temperature'] < 16 ]:
     ...     print(record)
-    (20.0, 15.0, 35.0, '2')
+    (np.float32(20.0), np.float32(15.0), np.float32(35.0), '2')
 
 But this is discouraged as this will be deprecated soon. The ``.iterdata()`` is
 therefore highly recommended.
@@ -633,19 +633,19 @@ class SequenceType(StructureType):
 
         >>> for line in seq.iterdata():
         ...     print(line)
-        (10, 15.2, 'Diamond_St')
-        (11, 13.1, 'Blacktail_Loop')
-        (12, 13.3, 'Platinum_St')
-        (13, 12.1, 'Kodiak_Trail')
+        (np.int32(10), np.float32(15.2), 'Diamond_St')
+        (np.int32(11), np.float32(13.1), 'Blacktail_Loop')
+        (np.int32(12), np.float32(13.3), 'Platinum_St')
+        (np.int32(13), np.float32(12.1), 'Kodiak_Trail')
 
     The order of the variables can be changed:
 
         >>> for line in seq['temperature', 'site', 'index'].iterdata():
         ...     print(line)
-        (15.2, 'Diamond_St', 10)
-        (13.1, 'Blacktail_Loop', 11)
-        (13.3, 'Platinum_St', 12)
-        (12.1, 'Kodiak_Trail', 13)
+        (np.float32(15.2), 'Diamond_St', np.int32(10))
+        (np.float32(13.1), 'Blacktail_Loop', np.int32(11))
+        (np.float32(13.3), 'Platinum_St', np.int32(12))
+        (np.float32(12.1), 'Kodiak_Trail', np.int32(13))
 
     We can iterate over children:
 
@@ -660,9 +660,9 @@ class SequenceType(StructureType):
 
         >>> for line in seq[ seq.index > 10 ].iterdata():
         ...     print(line)
-        (11, 13.1, 'Blacktail_Loop')
-        (12, 13.3, 'Platinum_St')
-        (13, 12.1, 'Kodiak_Trail')
+        (np.int32(11), np.float32(13.1), 'Blacktail_Loop')
+        (np.int32(12), np.float32(13.3), 'Platinum_St')
+        (np.int32(13), np.float32(12.1), 'Kodiak_Trail')
 
         >>> for line in seq[ seq.index > 10 ]['site'].iterdata():
         ...     print(line)
@@ -673,16 +673,16 @@ class SequenceType(StructureType):
         >>> for line in (seq['site', 'temperature'][seq.index > 10]
         ...              .iterdata()):
         ...     print(line)
-        ('Blacktail_Loop', 13.1)
-        ('Platinum_St', 13.3)
-        ('Kodiak_Trail', 12.1)
+        ('Blacktail_Loop', np.float32(13.1))
+        ('Platinum_St', np.float32(13.3))
+        ('Kodiak_Trail', np.float32(12.1))
 
     Or slice it:
 
         >>> for line in seq[::2].iterdata():
         ...     print(line)
-        (10, 15.2, 'Diamond_St')
-        (12, 13.3, 'Platinum_St')
+        (np.int32(10), np.float32(15.2), 'Diamond_St')
+        (np.int32(12), np.float32(13.3), 'Platinum_St')
 
         >>> for line in seq[ seq.index > 10 ][::2]['site'].iterdata():
         ...     print(line)
diff --git a/src/pydap/responses/dods.py b/src/pydap/responses/dods.py
index a7b73dc..1fc2d23 100644
--- a/src/pydap/responses/dods.py
+++ b/src/pydap/responses/dods.py
@@ -43,7 +43,9 @@ def DAP2_response_dtypemap(dtype):
 
 
 def tostring_with_byteorder(x, dtype):
-    return x.astype(dtype.str).newbyteorder(dtype.byteorder).tobytes()
+    order = dtype.byteorder
+    _var = x.astype(dtype.str)
+    return _var.view(_var.dtype.newbyteorder(order)).tobytes()
 
 
 class DODSResponse(BaseResponse):


Test patch:
diff --git a/.github/workflows/run_tests_macos.yml b/.github/workflows/run_tests_macos.yml
index 249a4bb..483c3de 100644
--- a/.github/workflows/run_tests_macos.yml
+++ b/.github/workflows/run_tests_macos.yml
@@ -27,7 +27,7 @@ jobs:
       run: |
         pip3 install -U pip setuptools
         pip3 install -U "werkzeug<2.1"
-        pip3 install -U "numpy<2.0" webtest gsw coards
+        pip3 install -U numpy webtest gsw coards
         pip3 install -U netCDF4 pytest Webob Jinja2 docopt-ng beautifulsoup4 requests testresources
     - name: Run tests with pytest
       run: |
diff --git a/.github/workflows/run_tests_ubuntu.yml b/.github/workflows/run_tests_ubuntu.yml
index 62ac12e..03aee57 100644
--- a/.github/workflows/run_tests_ubuntu.yml
+++ b/.github/workflows/run_tests_ubuntu.yml
@@ -25,7 +25,7 @@ jobs:
         python -m pip install --upgrade pip
         pip3 install -U pip setuptools
         pip3 install -U "werkzeug>=2.2.2" cryptography
-        pip3 install -U "numpy<2.0" webtest gsw coards
+        pip3 install -U numpy webtest gsw coards
         pip3 install -U netCDF4 pytest Webob Jinja2 docopt-ng beautifulsoup4 requests testresources
     - name: Run tests with pytest
       run: |
diff --git a/src/pydap/tests/test_responses_dods.py b/src/pydap/tests/test_responses_dods.py
index 2f72bbe..78c7145 100644
--- a/src/pydap/tests/test_responses_dods.py
+++ b/src/pydap/tests/test_responses_dods.py
@@ -5,6 +5,7 @@ import unittest
 from collections import OrderedDict
 
 import numpy as np
+from numpy.lib import Arrayterator
 from webob.headers import ResponseHeaders
 from webtest import TestApp as App
 
@@ -170,6 +171,14 @@ class TestDODSResponseStructure(unittest.TestCase):
 
     def test_body(self):
         """Test response body."""
+        arr = np.array(-10, dtype=np.int16).astype(">i")
+        arr_1d = np.array([arr])
+        newbyteo = arr_1d.dtype.newbyteorder(">")
+        newarr = arr_1d.astype(newbyteo)
+        # View with new byte order
+        expected_result = newarr.view(newbyteo)
+
+        # begins here
         dds, xdrdata = self.res.body.split(b"\nData:\n", 1)
         dds = dds.decode("ascii")
         self.assertEqual(
@@ -200,15 +209,7 @@ class TestDODSResponseStructure(unittest.TestCase):
                 ("i32", (np.array(-10, dtype=np.int32).astype(">i").tobytes())),
                 ("ui32", (np.array(10, dtype=np.uint32).astype(">I").tobytes())),
                 # -10 int16 upconv. to int32 for transfer
-                (
-                    "i16",
-                    (
-                        np.array(-10, dtype=np.int16)
-                        .astype(">i")
-                        .newbyteorder(">")
-                        .tobytes()
-                    ),
-                ),
+                ("i16", expected_result.tobytes()),
                 # 10 uint16 upconv. to uint32 for transfer
                 ("ui16", (np.array(10, dtype=np.uint16).astype(">I").tobytes())),
                 ("f32", (np.array(100, dtype=np.float32).astype(">f").tobytes())),
@@ -315,7 +316,7 @@ class TestDODSResponseArrayterator(unittest.TestCase):
     def test_grid(self):
         """Test a grid."""
         modified = copy.copy(SimpleGrid)
-        modified.SimpleGrid.SimpleGrid.data = np.lib.arrayterator.Arrayterator(
+        modified.SimpleGrid.SimpleGrid.data = Arrayterator(
             modified.SimpleGrid.SimpleGrid.data, 2
         )