Repository: ContactEngineering/SurfaceTopography (language: python)

Issue:
Error When Opening DATX Files: Unpacking Values Issue
## Description
When attempting to open `.datx` files (e.g. [test.datx](https://www.dropbox.com/scl/fi/zbhvolt2rip5qs4bbqiaf/test.datx?rlkey=u4k6wn4ybaaxuoulw423rpqzm&dl=0)) using the `open_topography` function, an error occurs related to unpacking values. This issue prevents successful file operations on `.datx` format data.

## Error Message
The following traceback details the specific error encountered:
```
Traceback (most recent call last):
  File "/Users/shirose/Library/CloudStorage/Dropbox/test.py", line 11, in <module>
    reader = open_topography(filename, format='datx')
  File "/Users/shirose/Library/Python/3.12/lib/python/site-packages/SurfaceTopography/IO/__init__.py", line 240, in open_topography
    return lookup_reader_by_format[format](fobj, **kwargs)
  File "/Users/shirose/Library/Python/3.12/lib/python/site-packages/SurfaceTopography/IO/DATX.py", line 101, in __init__
    category, unit, values = x_converter
ValueError: too many values to unpack (expected 3)
```


Discussion / hints:
I'm not an expert in DATX format or HDF5 technology. However, I wanted to share a Python code snippet, adapted by ChatGPT from a MATLAB sample provided by Zygo, that appears to successfully read DATX files. I hope this can help address the issue reported above.

```
import numpy as np
import matplotlib.pyplot as plt
import h5py

def get_z_data_with_units(map_hdf):
    """Convert Z data based on the unit converter attributes."""
    unit_out = 'um'  # Default output unit
    z_converter = map_hdf.attrs['Z Converter'][0]  # Correct access to the tuple inside the list
    category = z_converter[0].decode('utf-8')
    base_unit = z_converter[1].decode('utf-8')  # Correct decoding of base unit
    values = z_converter[2]  # Values array

    # Select the scaling based on the unit
    if base_unit == 'NanoMeters':
        z_scale = 1e-3
    elif base_unit == 'MicroMeters':
        z_scale = 1.0
    elif base_unit == 'MilliMeters':
        z_scale = 1e3
    elif base_unit == 'Meters':
        z_scale = 1e6
    elif base_unit == 'Fringes':
        S = values[2]
        O = values[3]
        W = values[1]
        z_scale = 1e6 * S * O * W
    else:
        unit_out = base_unit
        z_scale = 1

    print(base_unit)
    height_um = map_hdf[:] * z_scale  # Correct division using the scale
    return height_um, unit_out

def read_datx(file_name):
    with h5py.File(file_name, 'r') as file:
        # Find the first dataset under 'Data/Surface'
        base_path = 'Data/Surface'
        if base_path in file:
            dataset_names = list(file[base_path].keys())
            if dataset_names:
                dataset_path = f'{base_path}/{dataset_names[0]}'
                map_hdf = file[dataset_path]
                surface_data, unit_out = get_z_data_with_units(map_hdf)

                # Set visualization range
                vmin = -20e-3
                vmax = 20e-3

                # Visualization
                plt.figure()
                plt.imshow(surface_data, cmap='viridis', interpolation='none', vmin=vmin, vmax=vmax)
                plt.colorbar(label=f'Height [{unit_out}]')
                plt.title('Surface Data Visualization')
                plt.show()
            else:
                print("No datasets found under the specified base path.")
        else:
            print("Specified base path does not exist in the file.")

# Example usage
file_name = './test.datx'
read_datx(file_name)
```

Fix patch:
diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml
index 147aa1912..0c257cee9 100644
--- a/.github/workflows/wheels.yml
+++ b/.github/workflows/wheels.yml
@@ -18,7 +18,8 @@ jobs:
       matrix:
         buildplat:
         - [ubuntu-22.04, manylinux, x86_64]
-        # - [macos-10.15, macosx, x86_64]
+        - [ubuntu-22.04, manylinux, aarch64]
+        - [macos-14, macosx, arm64]
         # - [windows-2019, win, AMD64]
         python: ["cp38", "cp39", "cp310", "cp311", "cp312"]
       fail-fast: false
@@ -27,6 +28,11 @@ jobs:
       IS_32_BIT: ${{ matrix.buildplat[2] == 'x86' }}
 
     steps:
+      - if: matrix.buildplat[2] == 'aarch64'
+        uses: docker/setup-qemu-action@v3
+        with:
+          platforms: all
+
       - uses: actions/checkout@v4
         with:
           fetch-depth: 0
@@ -39,11 +45,11 @@ jobs:
       - name: Build wheels
         uses: pypa/cibuildwheel@v2.19.2
         env:
-          CIBW_BEFORE_ALL_LINUX: yum install -y hdf5 hdf5-devel netcdf netcdf-devel
           CIBW_BUILD: ${{ matrix.python }}-${{ matrix.buildplat[1] }}*
           CIBW_ARCHS: ${{ matrix.buildplat[2] }}
           CIBW_ENVIRONMENT_PASS_LINUX: RUNNER_OS
           CIBW_TEST_COMMAND: python -c 'import SurfaceTopography; print(SurfaceTopography.__version__)'
+          CIBW_TEST_SKIP: "*-*linux_{aarch64,ppc64le,s390x}"
 
       - uses: actions/upload-artifact@v4
         with:
diff --git a/SurfaceTopography/ChangeLog.md b/SurfaceTopography/ChangeLog.md
index f0996d3b9..e059d54bb 100644
--- a/SurfaceTopography/ChangeLog.md
+++ b/SurfaceTopography/ChangeLog.md
@@ -1,10 +1,19 @@
 Change log for SurfaceTopography
 ================================
 
+v1.15.0 (not yet released)
+--------------------------
+
+- ENH: Variable bandwidth analysis for topographies with missing data
+- ENH: NMM (Nanomeasuring Machine) reader
+- MAINT: Support for numpy >= 2.0.0
+- MAINT: Version discovery through `DiscoverVersion` module
+- CI: macOS and Linux aarch wheels
+
 v1.14.2 (24Jul24)
 -----------------
 
-- BUG: Fixed AL3D reader
+- BUG: Fixed data layout in AL3D reader
 
 v1.14.1 (16June24)
 ------------------
@@ -19,7 +28,7 @@ v1.14.0 (11June24)
 v1.13.13 (16Apr24)
 ------------------
 
-- BUG: Fixed reader for MetroPro v3 files, close #364
+- BUG: Fixed reader for MetroPro v3 files, closes #364
 - BUG: Wrong data layout in SUR files
 - TST: scipy.signal.triang -> scipy.signal.windows.triang, closes #365
 
diff --git a/SurfaceTopography/Container/IO/CE.py b/SurfaceTopography/Container/IO/CE.py
index 762460499..8dd0182e1 100644
--- a/SurfaceTopography/Container/IO/CE.py
+++ b/SurfaceTopography/Container/IO/CE.py
@@ -31,8 +31,8 @@
 import numpy as np
 import yaml
 
-from ...DiscoverVersion import __version__
 from ...IO import open_topography
+from ...Version import __version__
 from ..SurfaceContainer import LazySurfaceContainer, SurfaceContainer
 from .Reader import ContainerReaderBase
 
diff --git a/SurfaceTopography/DiscoverVersion.py b/SurfaceTopography/DiscoverVersion.py
deleted file mode 100644
index d9c3e29da..000000000
--- a/SurfaceTopography/DiscoverVersion.py
+++ /dev/null
@@ -1,84 +0,0 @@
-#
-# Copyright 2020-2023 Lars Pastewka
-#
-# ### MIT license
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-# SOFTWARE.
-#
-
-import os
-import subprocess
-
-
-class CannotDiscoverVersion(Exception):
-    pass
-
-
-def get_version_from_git():
-    """
-    Discover version from git repository.
-    """
-    if not os.path.exists('.git'):
-        raise CannotDiscoverVersion('.git subdirectory does not exist.')
-
-    try:
-        git_describe = subprocess.run(
-            ['git', 'describe', '--tags', '--dirty', '--always'],
-            stdout=subprocess.PIPE)
-    except FileNotFoundError:
-        git_describe = None
-    if git_describe is None or git_describe.returncode != 0:
-        raise CannotDiscoverVersion('git execution failed.')
-    version = git_describe.stdout.decode('latin-1').strip()
-
-    dirty = version.endswith('-dirty')
-
-    # Make version PEP 440 compliant
-    if dirty:
-        version = version.replace('-dirty', '')
-    version = version.strip('v')  # Remove leading 'v' if it exists
-    version = version.replace('-', '.dev', 1)
-    version = version.replace('-', '+', 1)
-    if dirty:
-        version += '.dirty'
-
-    return version
-
-
-_pkg_name = __name__.replace('.DiscoverVersion', '')
-
-# importlib is present in Python >= 3.8
-try:
-    from importlib.metadata import version
-
-    __version__ = version(_pkg_name)
-except ImportError:
-    __version__ = None
-
-# git works if we are in the source repository
-if __version__ is None:
-    try:
-        __version__ = get_version_from_git()
-    except CannotDiscoverVersion:
-        __version__ = None
-
-# Nope. Out of options.
-
-if __version__ is None:
-    raise CannotDiscoverVersion('Tried importlib and git')
diff --git a/SurfaceTopography/IO/DATX.py b/SurfaceTopography/IO/DATX.py
index 65245dc66..e4f8ad21b 100644
--- a/SurfaceTopography/IO/DATX.py
+++ b/SurfaceTopography/IO/DATX.py
@@ -98,19 +98,25 @@ def recursive_replace(raw_metadata, context):
                 z_converter, = h5[self._surface_path].attrs['Z Converter']
 
                 # Interpret converters and throw error if unsupported
-                category, unit, values = x_converter
+                category = x_converter[0]
+                unit = x_converter[1]
+                values = x_converter[-1]
                 if category != b'LateralCat' and unit != b'Pixels':
                     raise UnsupportedFormatFeature('DATX reader only supports `LateralCat` with `Pixels` unit for '
                                                    'X converter.')
                 physical_sizes_x = self._nb_grid_pts[0] * values[1]  # in units of meters!
 
-                category, unit, values = y_converter
+                category = y_converter[0]
+                unit = y_converter[1]
+                values = y_converter[-1]
                 if category != b'LateralCat' and unit != b'Pixels':
                     raise UnsupportedFormatFeature('DATX reader only supports `LateralCat` with `Pixels` unit for '
                                                    'Y converter.')
                 physical_sizes_y = self._nb_grid_pts[1] * values[1]  # in units of meters!
 
-                category, height_unit, values = z_converter
+                category = z_converter[0]
+                height_unit = z_converter[1]
+                values = z_converter[-1]
                 if category != b'HeightCat':
                     raise UnsupportedFormatFeature('DATX reader only supports `LateralCat` for Z converter.')
                 height_unit = mangle_length_unit_utf8(height_unit.decode('ascii'))
diff --git a/SurfaceTopography/IO/NC.py b/SurfaceTopography/IO/NC.py
index 61b810536..5cac8b7b9 100644
--- a/SurfaceTopography/IO/NC.py
+++ b/SurfaceTopography/IO/NC.py
@@ -413,6 +413,7 @@ def write_nc_uniform(topography, fobj, format='NETCDF3_64BIT_OFFSET'):
         mask_var = None
         if topography.has_undefined_data:
             mask_var = nc.createVariable('mask', 'b', spatial, **var_kwargs)
+            mask_var.maskandscale = False  # numpy-2: _FillValue overflows
 
         # Create variables for x- and y-positions, but only if physical_sizes
         # exist. (physical_sizes should always exist, but who knows...)
diff --git a/SurfaceTopography/IO/OPDx.py b/SurfaceTopography/IO/OPDx.py
index ecc849f26..34b30beb9 100644
--- a/SurfaceTopography/IO/OPDx.py
+++ b/SurfaceTopography/IO/OPDx.py
@@ -33,11 +33,13 @@
 import dateutil.parser
 import numpy as np
 
-from ..Exceptions import FileFormatMismatch, MetadataAlreadyFixedByFile, UnsupportedFormatFeature
+from ..Exceptions import (FileFormatMismatch, MetadataAlreadyFixedByFile,
+                          UnsupportedFormatFeature)
+from ..Support.UnitConversion import (get_unit_conversion_factor,
+                                      mangle_length_unit_utf8)
 from ..UniformLineScanAndTopography import Topography, UniformLineScan
-from ..Support.UnitConversion import get_unit_conversion_factor, mangle_length_unit_utf8
 from .common import OpenFromAny
-from .Reader import ReaderBase, ChannelInfo
+from .Reader import ChannelInfo, ReaderBase
 
 MAGIC = b'VCA DATA\x01\x00\x00\x55'
 
@@ -438,13 +440,13 @@ def read_dimension2d_content(stream):
 def _read_name(stream):
     # Names always have a size of 4 bytes
     length = _read_scalar(stream, '<u4')
-    return stream.read(length).decode('UTF-8')
+    return stream.read(length).decode('utf-8')
 
 
 def _read_string(stream):
     # String have variable lengths
     string_length = _read_varlen(stream)
-    return stream.read(string_length).decode('UTF-8')
+    return stream.read(string_length).decode('utf-8')
 
 
 def _read_scalar(stream, dtype):
@@ -459,13 +461,15 @@ def _read_scalar(stream, dtype):
 def _read_varlen(stream):
     lenlen = _read_scalar(stream, 'B')
     if lenlen == 1:
-        return _read_scalar(stream, 'B')
+        L = _read_scalar(stream, 'B')
     elif lenlen == 2:
-        return _read_scalar(stream, '<u2')
+        L = _read_scalar(stream, '<u2')
     elif lenlen == 4:
-        return _read_scalar(stream, '<u4')
+        L = _read_scalar(stream, '<u4')
     else:
         raise ValueError(f"Don't know how to read a variable length of size {lenlen}.")
+    # numpy-2: Cast to int because otherwise arithmetics will overflow
+    return int(L)
 
 
 def _read_string_list(stream):
diff --git a/SurfaceTopography/IO/XYZ.py b/SurfaceTopography/IO/XYZ.py
index 5d2b207e7..896ce6f71 100644
--- a/SurfaceTopography/IO/XYZ.py
+++ b/SurfaceTopography/IO/XYZ.py
@@ -409,7 +409,7 @@ def _read(self, physical_sizes=None, height_scale_factor=None, unit=None, info={
             x, z = data
         elif len(data) == 3:
             x, y, z = data
-            if ((x[1:] - x[0]) / (y[1:] - y[0])).ptp() < tol:
+            if np.ptp((x[1:] - x[0]) / (y[1:] - y[0])) < tol:
                 # This is a line scan and the first column is likely an index column.
                 kind = 'xz'
                 x = y
diff --git a/SurfaceTopography/Models/SelfAffine.py b/SurfaceTopography/Models/SelfAffine.py
index e3b49ae0c..d5f9a4f44 100644
--- a/SurfaceTopography/Models/SelfAffine.py
+++ b/SurfaceTopography/Models/SelfAffine.py
@@ -23,14 +23,15 @@
 # SOFTWARE.
 #
 import numpy as np
-from SurfaceTopography.Generation import fourier_synthesis
 
+from SurfaceTopography.Generation import fourier_synthesis
 
 # TODO abstract class for discrete and model PSDs
 # TODO all wavelengths and prefactors as properties
 #
 # TODO: think about abstract Statistical roughness hierarchy
 
+
 class AbstractStatisticalRoughness(object):
     def power_spectrum(seld):
         raise NotImplementedError
@@ -52,7 +53,7 @@ def __init__(self,
                  cr,
                  rolloff_wavelength,
                  hurst_exponent,
-                 longcut_wavelength=np.infty,
+                 longcut_wavelength=np.inf,
                  shortcut_wavelength=0, unit=None):
         self.cr = cr  # Ampliture of the PSD at and below the rolloff wavevector
         self.rolloff_wavelength = rolloff_wavelength
@@ -61,7 +62,7 @@ def __init__(self,
         self.hurst_exponent = hurst_exponent
 
         if self.shortcut_wavelength == 0:
-            self.shortcut_wavevector = np.infty
+            self.shortcut_wavevector = np.inf
         else:
             self.shortcut_wavevector = 2 * np.pi / self.shortcut_wavelength
 
@@ -129,7 +130,7 @@ def rms_height(self, shortcut_wavelength=None, longcut_wavelength=None, ):
         if longcut_wavelength is None:
             longcut_wavelength = self.longcut_wavelength
 
-        shortcut_wavevector = (2 * np.pi) / shortcut_wavelength if shortcut_wavelength > 0 else np.infty
+        shortcut_wavevector = (2 * np.pi) / shortcut_wavelength if shortcut_wavelength > 0 else np.inf
         longcut_wavevector = (2 * np.pi) / longcut_wavelength
 
         # see labbook of 220304
@@ -196,7 +197,7 @@ def variance_derivative(self, order, shortcut_wavelength=None, longcut_wavelengt
         else:
             longcut_wavelength = min(longcut_wavelength, self.longcut_wavelength)
 
-        shortcut_wavevector = (2 * np.pi) / shortcut_wavelength if shortcut_wavelength > 0 else np.infty
+        shortcut_wavevector = (2 * np.pi) / shortcut_wavelength if shortcut_wavelength > 0 else np.inf
         longcut_wavevector = (2 * np.pi) / longcut_wavelength
         rolloff_wavevector = self.rolloff_wavevector
         # prefactor of the self-afine region:
diff --git a/SurfaceTopography/Uniform/Detrending.py b/SurfaceTopography/Uniform/Detrending.py
index abec632cd..e2df268c6 100644
--- a/SurfaceTopography/Uniform/Detrending.py
+++ b/SurfaceTopography/Uniform/Detrending.py
@@ -140,7 +140,7 @@ def polyfit2_topography(self, full_output=False):
     x_grids = np.meshgrid(*x_grids, indexing='ij')
     # Quadratic terms
     x, y = x_grids
-    x_grids += [x * x, y * y, x * y]
+    x_grids += (x * x, y * y, x * y)
     if np.ma.getmask(arr) is np.ma.nomask:
         columns = [x.reshape((-1, 1)) for x in x_grids]
     else:
diff --git a/SurfaceTopography/Uniform/VariableBandwidth.py b/SurfaceTopography/Uniform/VariableBandwidth.py
index 0e6d06a21..151c710d5 100644
--- a/SurfaceTopography/Uniform/VariableBandwidth.py
+++ b/SurfaceTopography/Uniform/VariableBandwidth.py
@@ -87,7 +87,7 @@ def checkerboard_detrend_profile(self, subdivisions, order=1, return_plane=False
 
     b = np.array([np.bincount(region_index, h * (x ** i)) for i in range(order + 1)])
     C = np.array([[np.bincount(region_index, x ** (k + i)) for i in range(order + 1)] for k in range(order + 1)])
-    a = np.linalg.solve(C.T, b.T).T
+    a = np.linalg.solve(C.T, b.T.reshape(b.T.shape + (1,))).T[0]
 
     detrended_h = h - np.sum([a[i, region_index] * x ** i for i in range(order + 1)], axis=0)
     detrended_h.shape = shape
@@ -168,7 +168,7 @@ def checkerboard_detrend_area(self, subdivisions, order=1, return_plane=False):
 
     b = np.array([np.bincount(region_index, h * (x ** i) * (y ** j)) for i, j in ij])
     C = np.array([[np.bincount(region_index, x ** (k + i) * y ** (l + j)) for i, j in ij] for k, l in ij])
-    a = np.linalg.solve(C.T, b.T).T
+    a = np.linalg.solve(C.T, b.T.reshape(b.T.shape + (1,))).T[0]
 
     detrended_h = h - np.sum([a[k, region_index] * (x ** i) * (y ** j) for k, (i, j) in enumerate(ij)], axis=0)
     detrended_h.shape = shape
diff --git a/SurfaceTopography/Uniform/common.py b/SurfaceTopography/Uniform/common.py
index 43c6d6949..4499ce14e 100644
--- a/SurfaceTopography/Uniform/common.py
+++ b/SurfaceTopography/Uniform/common.py
@@ -32,8 +32,8 @@
 
 from ..FFTTricks import make_fft
 from ..HeightContainer import UniformTopographyInterface
-from ..UniformLineScanAndTopography import Topography
-from ..UniformLineScanAndTopography import DecoratedUniformTopography
+from ..UniformLineScanAndTopography import (DecoratedUniformTopography,
+                                            Topography)
 
 
 def bandwidth(self):
@@ -133,7 +133,7 @@ def plot(topography, subplot_location=111):
 
 
 class FilledTopography(DecoratedUniformTopography):
-    def __init__(self, topography, fill_value=-np.infty, info={}):
+    def __init__(self, topography, fill_value=-np.inf, info={}):
         """
         masked (undefined) data is replaced with `fill_value`.
 
diff --git a/SurfaceTopography/Version.py b/SurfaceTopography/Version.py
new file mode 100644
index 000000000..605709bc2
--- /dev/null
+++ b/SurfaceTopography/Version.py
@@ -0,0 +1,27 @@
+#
+# Copyright 2024 Lars Pastewka
+#
+# ### MIT license
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+#
+
+from DiscoverVersion import get_version
+
+__version__ = get_version('SurfaceTopography', __file__)
diff --git a/SurfaceTopography/__init__.py b/SurfaceTopography/__init__.py
index 8e4c3ca6c..1cb7fda6f 100644
--- a/SurfaceTopography/__init__.py
+++ b/SurfaceTopography/__init__.py
@@ -65,14 +65,12 @@
 
 from .Container import (open_container, read_container,  # noqa: F401
                         read_published_container)
-# This needs to be the first import, because __version__ is needed in the
-# packages that are subsequently imported.
-from .DiscoverVersion import __version__  # noqa: F401
 from .IO import open_topography, read_topography  # noqa: F401
 from .NonuniformLineScan import NonuniformLineScan  # noqa: F401
 from .Special import PlasticTopography, make_sphere  # noqa: F401
 from .UniformLineScanAndTopography import Topography  # noqa: F401
 from .UniformLineScanAndTopography import UniformLineScan  # noqa: F401
+from .Version import __version__  # noqa: F401
 
 # Add contact.engineering paper to bibliography
 SurfaceTopography.Support.Bibliography._default_dois = set(['10.1088/2051-672X/ac860a'])
diff --git a/SurfaceTopography/meson.build b/SurfaceTopography/meson.build
index cad42dff2..1a4ad7ae6 100644
--- a/SurfaceTopography/meson.build
+++ b/SurfaceTopography/meson.build
@@ -3,7 +3,7 @@
 # Pure Python sources
 python_sources = [
     '__init__.py',
-    'DiscoverVersion.py',
+    'Version.py',
     'Exceptions.py',
     'FFTTricks.py',
     'Generation.py',
diff --git a/c/autocorrelation.cpp b/c/autocorrelation.cpp
index a89be848c..cdab7422a 100644
--- a/c/autocorrelation.cpp
+++ b/c/autocorrelation.cpp
@@ -34,6 +34,7 @@ SOFTWARE.
 #include <cmath>
 
 #include <Python.h>
+#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION
 #define PY_ARRAY_UNIQUE_SYMBOL PYCO_ARRAY_API
 #define NO_IMPORT_ARRAY
 #include <numpy/arrayobject.h>
@@ -63,46 +64,46 @@ nonuniform_autocorrelation(PyObject *self, PyObject *args)
     FINALIZE_AND_RETURN;
   }
 
-  py_double_x = (PyObject*) PyArray_FROMANY((PyObject *) py_x, NPY_DOUBLE, 1, 1, NPY_C_CONTIGUOUS);
+  py_double_x = (PyObject*) PyArray_FROMANY((PyObject *) py_x, NPY_DOUBLE, 1, 1, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_double_x) {
     FINALIZE_AND_RETURN;
   }
 
-  py_double_h = (PyObject*) PyArray_FROMANY((PyObject *) py_h, NPY_DOUBLE, 1, 1, NPY_C_CONTIGUOUS);
+  py_double_h = (PyObject*) PyArray_FROMANY((PyObject *) py_h, NPY_DOUBLE, 1, 1, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_double_h) {
     FINALIZE_AND_RETURN;
   }
 
-  npy_intp nb_grid_pts = PyArray_DIM(py_double_x, 0);
-  if (PyArray_DIM(py_double_h, 0) != nb_grid_pts) {
+  npy_intp nb_grid_pts = PyArray_DIM((PyArrayObject *) py_double_x, 0);
+  if (PyArray_DIM((PyArrayObject *) py_double_h, 0) != nb_grid_pts) {
     PyErr_SetString(PyExc_TypeError, "x- and y-arrays must contain identical number of data points.");
   }
 
-  double *x = (double *) PyArray_DATA(py_double_x);
-  double *h = (double *) PyArray_DATA(py_double_h);
+  double *x = (double *) PyArray_DATA((PyArrayObject *) py_double_x);
+  double *h = (double *) PyArray_DATA((PyArrayObject *) py_double_h);
 
   double *distances;
   if (py_distances && py_distances != Py_None) {
-    py_double_distances = (PyObject*) PyArray_FROMANY((PyObject *) py_distances, NPY_DOUBLE, 1, 1, NPY_C_CONTIGUOUS);
+    py_double_distances = (PyObject*) PyArray_FROMANY((PyObject *) py_distances, NPY_DOUBLE, 1, 1, NPY_ARRAY_C_CONTIGUOUS);
     if (!py_double_distances) {
       FINALIZE_AND_RETURN;
     }
-    distances = (double *) PyArray_DATA(py_double_distances);
+    distances = (double *) PyArray_DATA((PyArrayObject *) py_double_distances);
   }
   else {
     py_double_distances = PyArray_EMPTY(1, &nb_grid_pts, NPY_DOUBLE, 0);
     if (!py_double_distances) {
       FINALIZE_AND_RETURN;
     }
-    distances = (double *) PyArray_DATA(py_double_distances);
+    distances = (double *) PyArray_DATA((PyArrayObject *) py_double_distances);
     /* Create distance array */
     for (int i = 0; i < nb_grid_pts; ++i)  distances[i] = i*physical_size/nb_grid_pts;
   }
 
-  npy_intp nb_distance_pts = PyArray_DIM(py_double_distances, 0);
+  npy_intp nb_distance_pts = PyArray_DIM((PyArrayObject *) py_double_distances, 0);
 
   py_acf = PyArray_ZEROS(1, &nb_distance_pts, NPY_DOUBLE, 0);
-  double *acf = (double *) PyArray_DATA(py_acf);
+  double *acf = (double *) PyArray_DATA((PyArrayObject *) py_acf);
 
   for (int i = 0; i < nb_grid_pts-1; ++i) {
     double x1 = x[i];
diff --git a/c/bicubic.cpp b/c/bicubic.cpp
index ab320368d..d837f4dd4 100644
--- a/c/bicubic.cpp
+++ b/c/bicubic.cpp
@@ -35,6 +35,7 @@ SOFTWARE.
 #include <iostream>
 
 #include <Python.h>
+#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION
 #define PY_ARRAY_UNIQUE_SYMBOL PYCO_ARRAY_API
 #define NO_IMPORT_ARRAY
 #include <numpy/arrayobject.h>
@@ -357,14 +358,14 @@ bicubic_init(bicubic_t *self, PyObject *args, PyObject *kwargs)
   py_values = PyArray_FROMANY(py_values_in, NPY_DOUBLE, 2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_values)
     return -1;
-  nx = PyArray_DIM(py_values, 0);
-  ny = PyArray_DIM(py_values, 1);
+  nx = PyArray_DIM((PyArrayObject *) py_values, 0);
+  ny = PyArray_DIM((PyArrayObject *) py_values, 1);
 
   if (py_derivativex_in) {
     py_derivativex = PyArray_FROMANY(py_derivativex_in, NPY_DOUBLE, 2, 2, NPY_ARRAY_C_CONTIGUOUS);
     if (!py_derivativex)
       return -1;
-    if (PyArray_DIM(py_derivativex, 0) != nx || PyArray_DIM(py_derivativex, 1) != ny) {
+    if (PyArray_DIM((PyArrayObject *) py_derivativex, 0) != nx || PyArray_DIM((PyArrayObject *) py_derivativex, 1) != ny) {
 	  PyErr_SetString(PyExc_ValueError, "x-derivative must have same shape as values.");
 	  return -1;
     }
@@ -374,7 +375,7 @@ bicubic_init(bicubic_t *self, PyObject *args, PyObject *kwargs)
     py_derivativey = PyArray_FROMANY(py_derivativey_in, NPY_DOUBLE, 2, 2, NPY_ARRAY_C_CONTIGUOUS);
     if (!py_derivativey)
       return -1;
-    if (PyArray_DIM(py_derivativey, 0) != nx || PyArray_DIM(py_derivativey, 1) != ny) {
+    if (PyArray_DIM((PyArrayObject *) py_derivativey, 0) != nx || PyArray_DIM((PyArrayObject *) py_derivativey, 1) != ny) {
 	  PyErr_SetString(PyExc_ValueError, "y-derivative must have same shape as values.");
 	  return -1;
     }
@@ -382,18 +383,18 @@ bicubic_init(bicubic_t *self, PyObject *args, PyObject *kwargs)
 
   double *derivativex{NULL}, *derivativey{NULL};
   if (py_derivativex) {
-    derivativex = static_cast<double*>(PyArray_DATA(py_derivativex));
+    derivativex = static_cast<double*>(PyArray_DATA((PyArrayObject *) py_derivativex));
   }
   if (py_derivativey) {
-    derivativey = static_cast<double*>(PyArray_DATA(py_derivativey));
+    derivativey = static_cast<double*>(PyArray_DATA((PyArrayObject *) py_derivativey));
   }
 
   if (derivativex || derivativey) {
-    self->map_ = new Bicubic(nx, ny, static_cast<double*>(PyArray_DATA(py_values)), derivativex, derivativey, false,
+    self->map_ = new Bicubic(nx, ny, static_cast<double*>(PyArray_DATA((PyArrayObject *) py_values)), derivativex, derivativey, false,
                              false);
   }
   else {
-    self->map_ = new Bicubic(nx, ny, static_cast<double*>(PyArray_DATA(py_values)), NULL, NULL, true, false);
+    self->map_ = new Bicubic(nx, ny, static_cast<double*>(PyArray_DATA((PyArrayObject *) py_values)), NULL, NULL, true, false);
   }
   Py_DECREF(py_values);
   if (py_derivativex)  Py_DECREF(py_derivativex);
@@ -434,16 +435,16 @@ bicubic_call(bicubic_t *self, PyObject *args, PyObject *kwargs)
     if (!py_r)
       return NULL;
 
-    if (PyArray_DIM(py_r, 1) != 2) {
+    if (PyArray_DIM((PyArrayObject *) py_r, 1) != 2) {
       PyErr_SetString(PyExc_ValueError, "Map index needs to have x- and y-component only.");
       return NULL;
     }
 
-    npy_intp n = PyArray_DIM(py_r, 0);
-    double *r = (double *) PyArray_DATA(py_r);
+    npy_intp n = PyArray_DIM((PyArrayObject *) py_r, 0);
+    double *r = (double *) PyArray_DATA((PyArrayObject *) py_r);
 
     PyObject *py_v = PyArray_SimpleNew(1, &n, NPY_DOUBLE);
-    double *v = (double *) PyArray_DATA(py_v);
+    double *v = (double *) PyArray_DATA((PyArrayObject *) py_v);
 
     for (int i = 0; i < n; i++) {
       self->map_->eval(r[2*i], r[2*i+1], v[i]);
@@ -471,17 +472,17 @@ bicubic_call(bicubic_t *self, PyObject *args, PyObject *kwargs)
       return NULL;
 
     /* Check that x and y have the same number of dimensions */
-    if (PyArray_NDIM(py_xd) != PyArray_NDIM(py_yd)) {
+    if (PyArray_NDIM((PyArrayObject *) py_xd) != PyArray_NDIM((PyArrayObject *) py_yd)) {
       PyErr_SetString(PyExc_ValueError, "x- and y-components need to have identical number of dimensions.");
       return NULL;
     }
 
     /* Check that x and y have the same length in each dimension */
-    int ndims = PyArray_NDIM(py_xd);
-    npy_intp *dims = PyArray_DIMS(py_xd);
+    int ndims = PyArray_NDIM((PyArrayObject *) py_xd);
+    npy_intp *dims = PyArray_DIMS((PyArrayObject *) py_xd);
     npy_intp n = 1;
     for (int i = 0; i < ndims; i++) {
-      npy_intp d = PyArray_DIM(py_yd, i);
+      npy_intp d = PyArray_DIM((PyArrayObject *) py_yd, i);
 
       if (dims[i] != d) {
 	    PyErr_SetString(PyExc_ValueError, "x- and y-components vectors need to have the same length.");
@@ -491,26 +492,26 @@ bicubic_call(bicubic_t *self, PyObject *args, PyObject *kwargs)
       n *= d;
     }
 
-    double *x{static_cast<double *>(PyArray_DATA(py_xd))};
-    double *y{static_cast<double *>(PyArray_DATA(py_yd))};
+    double *x{static_cast<double *>(PyArray_DATA((PyArrayObject *) py_xd))};
+    double *y{static_cast<double *>(PyArray_DATA((PyArrayObject *) py_yd))};
 
     PyObject *py_v{PyArray_SimpleNew(ndims, dims, NPY_DOUBLE)};
-    double *v{static_cast<double *>(PyArray_DATA(py_v))};
+    double *v{static_cast<double *>(PyArray_DATA((PyArrayObject *) py_v))};
 
     PyObject *py_return_value;
 
     if (derivative > 0) {
       PyObject *py_dx{PyArray_SimpleNew(ndims, dims, NPY_DOUBLE)};
       PyObject *py_dy{PyArray_SimpleNew(ndims, dims, NPY_DOUBLE)};
-      double *dx{static_cast<double *>(PyArray_DATA(py_dx))};
-      double *dy{static_cast<double *>(PyArray_DATA(py_dy))};
+      double *dx{static_cast<double *>(PyArray_DATA((PyArrayObject *) py_dx))};
+      double *dy{static_cast<double *>(PyArray_DATA((PyArrayObject *) py_dy))};
       if (derivative > 1) {
         PyObject *py_d2x{PyArray_SimpleNew(ndims, dims, NPY_DOUBLE)};
         PyObject *py_d2y{PyArray_SimpleNew(ndims, dims, NPY_DOUBLE)};
         PyObject *py_d2xy{PyArray_SimpleNew(ndims, dims, NPY_DOUBLE)};
-        double *d2x{static_cast<double *>(PyArray_DATA(py_d2x))};
-        double *d2y{static_cast<double *>(PyArray_DATA(py_d2y))};
-        double *d2xy{static_cast<double *>(PyArray_DATA(py_d2xy))};
+        double *d2x{static_cast<double *>(PyArray_DATA((PyArrayObject *) py_d2x))};
+        double *d2y{static_cast<double *>(PyArray_DATA((PyArrayObject *) py_d2y))};
+        double *d2xy{static_cast<double *>(PyArray_DATA((PyArrayObject *) py_d2xy))};
         for (int i = 0; i < n; i++) {
           self->map_->eval(x[i], y[i], v[i], dx[i], dy[i], d2x[i], d2y[i], d2xy[i]);
         }
diff --git a/c/patchfinder.cpp b/c/patchfinder.cpp
index 1a1fe0f84..53002f9c0 100644
--- a/c/patchfinder.cpp
+++ b/c/patchfinder.cpp
@@ -34,6 +34,7 @@ SOFTWARE.
 #include <cstddef>
 
 #include <Python.h>
+#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION
 #define PY_ARRAY_UNIQUE_SYMBOL PYCO_ARRAY_API
 #define NO_IMPORT_ARRAY
 #include <numpy/arrayobject.h>
@@ -127,14 +128,14 @@ PyObject *assign_patch_numbers(PyObject *self, PyObject *args)
   if (py_stencil) {
     py_long_stencil =
       (PyArrayObject*) PyArray_FROMANY((PyObject *) py_stencil, NPY_LONG,
-                                       2, 2, NPY_C_CONTIGUOUS);
+                                       2, 2, NPY_ARRAY_C_CONTIGUOUS);
     if (!py_long_stencil)
       return NULL;
 
     sx = PyArray_DIM(py_long_stencil, 0);
     npy_intp sy = PyArray_DIM(py_long_stencil, 1);
 
-    stencil = (npy_long*) PyArray_DATA(py_long_stencil);
+    stencil = (npy_long*) PyArray_DATA((PyArrayObject *) py_long_stencil);
 
     if (sy != 2) {
       PyErr_SetString(PyExc_TypeError, "Stencil must have dimension 2 in the "
@@ -147,21 +148,21 @@ PyObject *assign_patch_numbers(PyObject *self, PyObject *args)
   }
 
   py_bool_map = (PyArrayObject*) PyArray_FROMANY((PyObject *) py_map, NPY_BOOL,
-                                                 2, 2, NPY_C_CONTIGUOUS);
+                                                 2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_bool_map)
     return NULL;
 
   npy_intp nx = PyArray_DIM(py_bool_map, 0);
   npy_intp ny = PyArray_DIM(py_bool_map, 1);
 
-  npy_bool *map = (npy_bool*) PyArray_DATA(py_bool_map);
+  npy_bool *map = (npy_bool*) PyArray_DATA((PyArrayObject *) py_bool_map);
 
   npy_intp dims[2] = { nx, ny };
 
   PyObject *py_id = PyArray_ZEROS(2, dims, NPY_INT, 0);
   if (!py_id)
     return NULL;
-  npy_int *id = (npy_int *) PyArray_DATA(py_id);
+  npy_int *id = (npy_int *) PyArray_DATA((PyArrayObject *) py_id);
 
   std::ptrdiff_t i, j;
   std::ptrdiff_t k = 0;
@@ -236,21 +237,21 @@ PyObject *assign_segment_numbers(PyObject *self, PyObject *args)
 
   PyArrayObject *py_bool_map = NULL;
   py_bool_map = (PyArrayObject*) PyArray_FROMANY((PyObject *) py_map, NPY_BOOL,
-                         2, 2, NPY_C_CONTIGUOUS);
+                         2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_bool_map)
     return NULL;
 
   npy_intp nx = PyArray_DIM(py_bool_map, 0);
   npy_intp ny = PyArray_DIM(py_bool_map, 1);
 
-  npy_bool *map = (npy_bool*) PyArray_DATA(py_bool_map);
+  npy_bool *map = (npy_bool*) PyArray_DATA((PyArrayObject *) py_bool_map);
 
   npy_intp dims[2] = { nx, ny };
 
   PyObject *py_id = PyArray_ZEROS(2, dims, NPY_INT, 0);
   if (!py_id)
     return NULL;
-  npy_int *id = (npy_int *) PyArray_DATA(py_id);
+  npy_int *id = (npy_int *) PyArray_DATA((PyArrayObject *) py_id);
 
   std::ptrdiff_t i, j, k = 0;
   npy_int p = 0;
@@ -288,15 +289,15 @@ PyObject *shortest_distance(PyObject *self, PyObject *args)
   PyArrayObject *py_bool_to = NULL;
 
   py_bool_fromc = (PyArrayObject*)
-    PyArray_FROMANY((PyObject *) py_fromc, NPY_BOOL, 2, 2, NPY_C_CONTIGUOUS);
+    PyArray_FROMANY((PyObject *) py_fromc, NPY_BOOL, 2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_bool_fromc)
     return NULL;
   py_bool_fromp = (PyArrayObject*)
-    PyArray_FROMANY((PyObject *) py_fromp, NPY_BOOL, 2, 2, NPY_C_CONTIGUOUS);
+    PyArray_FROMANY((PyObject *) py_fromp, NPY_BOOL, 2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_bool_fromp)
     return NULL;
   py_bool_to = (PyArrayObject*)
-    PyArray_FROMANY((PyObject *) py_to, NPY_BOOL, 2, 2, NPY_C_CONTIGUOUS);
+    PyArray_FROMANY((PyObject *) py_to, NPY_BOOL, 2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_bool_to)
     return NULL;
 
@@ -315,15 +316,15 @@ PyObject *shortest_distance(PyObject *self, PyObject *args)
     return NULL;
   }
 
-  npy_bool *fromc = (npy_bool*) PyArray_DATA(py_bool_fromc);
-  npy_bool *fromp = (npy_bool*) PyArray_DATA(py_bool_fromp);
-  npy_bool *to = (npy_bool*) PyArray_DATA(py_bool_to);
+  npy_bool *fromc = (npy_bool*) PyArray_DATA((PyArrayObject *) py_bool_fromc);
+  npy_bool *fromp = (npy_bool*) PyArray_DATA((PyArrayObject *) py_bool_fromp);
+  npy_bool *to = (npy_bool*) PyArray_DATA((PyArrayObject *) py_bool_to);
 
   npy_intp dims[2] = { nx, ny };
   PyObject *py_dist = PyArray_ZEROS(2, dims, NPY_DOUBLE, 0);
   if (!py_dist)
     return NULL;
-  npy_double *dist = (npy_double *) PyArray_DATA(py_dist);
+  npy_double *dist = (npy_double *) PyArray_DATA((PyArrayObject *) py_dist);
 
   /*
    * Make sure there is something to find
@@ -533,27 +534,27 @@ PyObject *distance_map(PyObject *self, PyObject *args)
   PyArrayObject *py_bool_map_xy = NULL;
 
   py_bool_map_xy = (PyArrayObject*)
-    PyArray_FROMANY((PyObject *) py_map_xy, NPY_BOOL, 2, 2, NPY_C_CONTIGUOUS);
+    PyArray_FROMANY((PyObject *) py_map_xy, NPY_BOOL, 2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_bool_map_xy)
     return NULL;
 
   npy_intp nx = PyArray_DIM(py_bool_map_xy, 0);
   npy_intp ny = PyArray_DIM(py_bool_map_xy, 1);
 
-  npy_bool *map_xy = (npy_bool*) PyArray_DATA(py_bool_map_xy);
+  npy_bool *map_xy = (npy_bool*) PyArray_DATA((PyArrayObject *) py_bool_map_xy);
 
   /* This stores the distance to the closest point on the contour */
   npy_intp dims[2] = { nx, ny };
   PyObject *py_dist_xy = PyArray_ZEROS(2, dims, NPY_DOUBLE, 0);
   if (!py_dist_xy)
     return NULL;
-  npy_double *dist_xy = (npy_double *) PyArray_DATA(py_dist_xy);
+  npy_double *dist_xy = (npy_double *) PyArray_DATA((PyArrayObject *) py_dist_xy);
 
   /* This stores the index of the closest point */
   PyObject *py_next_xy = PyArray_ZEROS(2, dims, NPY_INT, 0);
   if (!py_next_xy)
     return NULL;
-  npy_int *next_xy = (npy_int *) PyArray_DATA(py_next_xy);
+  npy_int *next_xy = (npy_int *) PyArray_DATA((PyArrayObject *) py_next_xy);
 
   /*
    * Fill map with maximum distance
@@ -676,27 +677,27 @@ PyObject *closest_patch_map(PyObject *self, PyObject *args)
   PyArrayObject *py_int_map_xy = NULL;
 
   py_int_map_xy = (PyArrayObject*)
-    PyArray_FROMANY((PyObject *) py_map_xy, NPY_INT, 2, 2, NPY_C_CONTIGUOUS);
+    PyArray_FROMANY((PyObject *) py_map_xy, NPY_INT, 2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_int_map_xy)
     return NULL;
 
   npy_intp nx = PyArray_DIM(py_int_map_xy, 0);
   npy_intp ny = PyArray_DIM(py_int_map_xy, 1);
 
-  npy_int *map_xy = (npy_int*) PyArray_DATA(py_int_map_xy);
+  npy_int *map_xy = (npy_int*) PyArray_DATA((PyArrayObject *) py_int_map_xy);
 
   /* This stores the distance to the closest point on the contour */
   npy_intp dims[2] = { nx, ny };
   PyObject *py_dist_xy = PyArray_ZEROS(2, dims, NPY_DOUBLE, 0);
   if (!py_dist_xy)
     return NULL;
-  npy_double *dist_xy = (npy_double *) PyArray_DATA(py_dist_xy);
+  npy_double *dist_xy = (npy_double *) PyArray_DATA((PyArrayObject *) py_dist_xy);
 
   /* This stores the tag of the closest point */
   PyObject *py_next_xy = PyArray_ZEROS(2, dims, NPY_INT, 0);
   if (!py_next_xy)
     return NULL;
-  npy_int *next_xy = (npy_int *) PyArray_DATA(py_next_xy);
+  npy_int *next_xy = (npy_int *) PyArray_DATA((PyArrayObject *) py_next_xy);
 
   /*
    * Fill map with maximum distance
@@ -800,11 +801,11 @@ PyObject *correlation_function(PyObject *self, PyObject *args)
   PyArrayObject *py_double_map1 = NULL, *py_double_map2 = NULL;
 
   py_double_map1 = (PyArrayObject*)
-    PyArray_FROMANY((PyObject *) py_map1, NPY_DOUBLE, 2, 2, NPY_C_CONTIGUOUS);
+    PyArray_FROMANY((PyObject *) py_map1, NPY_DOUBLE, 2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_double_map1)
     return NULL;
   py_double_map2 = (PyArrayObject*)
-    PyArray_FROMANY((PyObject *) py_map2, NPY_DOUBLE, 2, 2, NPY_C_CONTIGUOUS);
+    PyArray_FROMANY((PyObject *) py_map2, NPY_DOUBLE, 2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_double_map2)
     return NULL;
 
@@ -817,8 +818,8 @@ PyObject *correlation_function(PyObject *self, PyObject *args)
             "Both maps need to have the identical dimensions.");
   }
 
-  npy_double *map1 = (npy_double*) PyArray_DATA(py_double_map1);
-  npy_double *map2 = (npy_double*) PyArray_DATA(py_double_map2);
+  npy_double *map1 = (npy_double*) PyArray_DATA((PyArrayObject *) py_double_map1);
+  npy_double *map2 = (npy_double*) PyArray_DATA((PyArrayObject *) py_double_map2);
 
   /*
    * Correlation function
@@ -827,7 +828,7 @@ PyObject *correlation_function(PyObject *self, PyObject *args)
   PyObject *py_c = PyArray_ZEROS(1, dims, NPY_DOUBLE, 0);
   if (!py_c)
     return NULL;
-  npy_double *c = (npy_double*) PyArray_DATA(py_c);
+  npy_double *c = (npy_double*) PyArray_DATA((PyArrayObject *) py_c);
 
   /*
    * Number of points found at a certain distance
@@ -835,7 +836,7 @@ PyObject *correlation_function(PyObject *self, PyObject *args)
   PyObject *py_n = PyArray_ZEROS(1, dims, NPY_INT, 0);
   if (!py_n)
     return NULL;
-  npy_int *n = (npy_int*) PyArray_DATA(py_n);
+  npy_int *n = (npy_int*) PyArray_DATA((PyArrayObject *) py_n);
 
   /*
    * Fill with zeros
@@ -892,17 +893,17 @@ PyObject *correlation_function(PyObject *self, PyObject *args)
   PyObject *py_r = PyArray_ZEROS(1, dims, NPY_DOUBLE, 0);
   if (!py_r)
     return NULL;
-  npy_double *r = (npy_double*) PyArray_DATA(py_r);
+  npy_double *r = (npy_double*) PyArray_DATA((PyArrayObject *) py_r);
   /* Correlation function */
   PyObject *py_cc = PyArray_ZEROS(1, dims, NPY_DOUBLE, 0);
   if (!py_cc)
     return NULL;
-  npy_double *cc = (npy_double*) PyArray_DATA(py_cc);
+  npy_double *cc = (npy_double*) PyArray_DATA((PyArrayObject *) py_cc);
   /* Integrated correlation function */
   PyObject *py_Icc = PyArray_ZEROS(1, dims, NPY_DOUBLE, 0);
   if (!py_Icc)
     return NULL;
-  npy_double *Icc = (npy_double*) PyArray_DATA(py_Icc);
+  npy_double *Icc = (npy_double*) PyArray_DATA((PyArrayObject *) py_Icc);
 
   /*
    * Normalize and integrate
@@ -1012,14 +1013,14 @@ PyObject *perimeter_length(PyObject *self, PyObject *args)
   PyArrayObject *py_bool_map = NULL;
 
   py_bool_map = (PyArrayObject*) PyArray_FROMANY((PyObject *) py_map, NPY_BOOL,
-                         2, 2, NPY_C_CONTIGUOUS);
+                         2, 2, NPY_ARRAY_C_CONTIGUOUS);
   if (!py_bool_map)
     return NULL;
 
   npy_intp nx = PyArray_DIM(py_bool_map, 0);
   npy_intp ny = PyArray_DIM(py_bool_map, 1);
 
-  npy_bool *map = (npy_bool*) PyArray_DATA(py_bool_map);
+  npy_bool *map = (npy_bool*) PyArray_DATA((PyArrayObject *) py_bool_map);
 
   double length = 0.0;
 
diff --git a/discover_version.py b/discover_version.py
index 0f6045cee..a5cd03149 100644
--- a/discover_version.py
+++ b/discover_version.py
@@ -1,5 +1,5 @@
 #
-# Copyright 2022 Lars Pastewka
+# Copyright 2022-2024 Lars Pastewka
 #
 # ### MIT license
 #
@@ -22,73 +22,6 @@
 # SOFTWARE.
 #
 
-#
-# This is the most minimal-idiotic way of discovering the version that I
-# could come up with. It deals with the following issues:
-# * If we are installed, we can get the version from package metadata,
-#   either via importlib.metadata or from pkg_resources. This also holds for
-#   wheels that contain the metadata. We are good! Yay!
-# * If we are not installed, there are two options:
-#   - We are working within the source git repository. Then
-#        git describe --tags --always
-#     yields a reasonable version descriptor, but that is unfortunately not
-#     PEP 440 compliant (see https://peps.python.org/pep-0440/). We need to
-#     mangle the version string to yield something compatible.
-# - If we install from a source tarball, we need to parse PKG-INFO manually.
-#
-
-import subprocess
-
-
-class CannotDiscoverVersion(Exception):
-    pass
-
-
-def get_version_from_pkg_info():
-    """
-    Discover version from PKG-INFO file.
-    """
-    f = open('PKG-INFO', 'r')
-    l = f.readline()
-    while l:
-        if l.startswith('Version:'):
-            return l[8:].strip()
-        l = f.readline()
-    raise CannotDiscoverVersion("No line starting with 'Version:' in 'PKG-INFO'.")
-
-
-def get_version_from_git():
-    """
-    Discover version from git repository.
-    """
-    git_describe = subprocess.run(
-        ['git', 'describe', '--tags', '--dirty', '--always'],
-        stdout=subprocess.PIPE)
-    if git_describe.returncode != 0:
-        raise CannotDiscoverVersion('git execution failed.')
-    version = git_describe.stdout.decode('latin-1').strip()
-
-    dirty = version.endswith('-dirty')
-
-    # Make version PEP 440 compliant
-    if dirty:
-        version = version.replace('-dirty', '')
-    version = version.strip('v')  # Remove leading 'v' if it exists
-    version = version.replace('-', '.dev', 1)
-    version = version.replace('-', '+', 1)
-    if dirty:
-        version += '.dirty'
-
-    return version
-
-
-try:
-    version = get_version_from_git()
-except CannotDiscoverVersion:
-    version = get_version_from_pkg_info()
-
-#
-# Print version to screen
-#
+from DiscoverVersion import get_version
 
-print(version)
+print(get_version('SurfaceTopography', __file__))
diff --git a/pyproject.toml b/pyproject.toml
index 4598e36df..467580b55 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,5 +1,12 @@
 [build-system]
-requires = ["meson>=1.0.0", "meson-python>=0.16.0", "ninja", "numpy<2.0"]
+requires = [
+    "meson>=1.0.0",
+    "meson-python>=0.13.0",
+    "ninja",
+    "oldest-supported-numpy; python_version=='3.8'",
+    "numpy>=2.0.0; python_version>='3.9'",
+    "DiscoverVersion"
+]
 build-backend = "mesonpy"
 
 [project]
@@ -22,14 +29,15 @@ classifiers = [
 requires-python = ">=3.8.0"
 dynamic = [ "version" ]
 dependencies = [
-    "numpy<2.0",
+    "DiscoverVersion",
+    "numpy",
     "scipy>=1.9.0",
     "defusedxml",
     "h5py",
     "igor2",
     "matplotlib>=1.0.0",
     "muFFT>=0.91.2",
-    "NuMPI>=0.3.1",
+    "NuMPI>=0.6.1",
     "numpyencoder",
     "openpyxl",
     "Pillow",


Test patch:
diff --git a/.github/workflows/test-code-functionality.yml b/.github/workflows/test-code-functionality.yml
index 28c7fb615..98a9753fa 100644
--- a/.github/workflows/test-code-functionality.yml
+++ b/.github/workflows/test-code-functionality.yml
@@ -58,6 +58,8 @@ jobs:
             pip-compile pyproject.toml
             NUMPY_VERSION=$(grep 'numpy==' requirements.txt)        
             MUFFT_VERSION=$(grep 'mufft==' requirements.txt)
+            echo "numpy version: ${NUMPY_VERSION}"
+            echo "muFFT version: ${MUFFT_VERSION}"
             # Install numpy
             pip install -v ${NUMPY_VERSION}
             # Install mpi4py
@@ -84,7 +86,7 @@ jobs:
       - name: Install extension module
         run: |
           source venv/bin/activate
-          python3 -m pip install .[test]
+          python3 -m pip install --no-binary netCDF4 .[test]
 
       - name: Test with pytest
         run: |
diff --git a/test/IO/test_datx.py b/test/IO/test_datx.py
index 65550f8f4..f4c7f820f 100644
--- a/test/IO/test_datx.py
+++ b/test/IO/test_datx.py
@@ -26,7 +26,6 @@
 
 import numpy as np
 import pytest
-
 from NuMPI import MPI
 
 from SurfaceTopography.IO import DATXReader
@@ -36,7 +35,7 @@
     reason="tests only serial funcionalities, please execute with pytest")
 
 
-def test_datx_metadata(file_format_examples):
+def test_datx1_metadata(file_format_examples):
     file_path = os.path.join(file_format_examples, 'datx-1.datx')
 
     r = DATXReader(file_path)
@@ -56,3 +55,25 @@ def test_datx_metadata(file_format_examples):
     np.testing.assert_allclose(t.min(), -23393.847656, rtol=1e-6)
 
     np.testing.assert_allclose(t.rms_height_from_area(), 6304.986277, rtol=1e-6)
+
+
+def test_datx2_metadata(file_format_examples):
+    file_path = os.path.join(file_format_examples, 'datx-2.datx')
+
+    r = DATXReader(file_path)
+    t = r.topography()
+
+    nx, ny = t.nb_grid_pts
+    assert nx == 1000
+    assert ny == 1000
+
+    assert t.unit == 'nm'
+
+    sx, sy = t.physical_sizes
+    np.testing.assert_allclose(sx, 173108.716695, rtol=1e-6)
+    np.testing.assert_allclose(sy, 173108.716695, rtol=1e-6)
+
+    np.testing.assert_allclose(t.max(), 11.037771, rtol=1e-6)
+    np.testing.assert_allclose(t.min(), -116.920823, rtol=1e-6)
+
+    np.testing.assert_allclose(t.rms_height_from_area(), 1.905869, rtol=1e-6)
diff --git a/test/IO/test_io.py b/test/IO/test_io.py
index 036102605..22cd56e09 100644
--- a/test/IO/test_io.py
+++ b/test/IO/test_io.py
@@ -122,6 +122,7 @@ def _convert_filelist(filelist):
                                               'lext-1.lext',
                                               'lext-2.lext',
                                               'datx-1.datx',
+                                              'datx-2.datx',
                                               'oir-1.oir',
                                               'poir-1.poir',
                                               'stp-1.stp',
diff --git a/test/file_format_examples/datx-2.datx b/test/file_format_examples/datx-2.datx
new file mode 100644
index 000000000..f913dcb20
Binary files /dev/null and b/test/file_format_examples/datx-2.datx differ
diff --git a/test/test_parallel.py b/test/test_parallel.py
index c9e9f01bc..355f49bfd 100644
--- a/test/test_parallel.py
+++ b/test/test_parallel.py
@@ -27,7 +27,6 @@
 """Specific test for MPI-parallel functionality"""
 
 import numpy as np
-
 from NuMPI.Tools import Reduction
 
 from SurfaceTopography import Topography
@@ -69,7 +68,7 @@ def test_fill_undefined_data_parallel(comm):
                             subdomain_locations=(0, comm.rank)
                             )
 
-    filled_topography = topography.fill_undefined_data(fill_value=-np.infty)
+    filled_topography = topography.fill_undefined_data(fill_value=-np.inf)
     assert topography.has_undefined_data
     assert not filled_topography.has_undefined_data
 
@@ -79,5 +78,5 @@ def test_fill_undefined_data_parallel(comm):
     reduction = Reduction(comm)
 
     assert reduction.all(filled_topography[nmask] == topography[nmask])
-    assert reduction.all(filled_topography[mask] == - np.infty)
+    assert reduction.all(filled_topography[mask] == - np.inf)
     assert not filled_topography.has_undefined_data
diff --git a/test/test_pipeline.py b/test/test_pipeline.py
index c9ab3c1ac..7b3c4d42f 100644
--- a/test/test_pipeline.py
+++ b/test/test_pipeline.py
@@ -135,7 +135,7 @@ def test_fill_undefined_data_linescan(periodic):
                                  )
     assert topography.has_undefined_data
 
-    filled_topography = topography.fill_undefined_data(fill_value=-np.infty)
+    filled_topography = topography.fill_undefined_data(fill_value=-np.inf)
     assert not filled_topography.has_undefined_data
 
     assert filled_topography.physical_sizes == topography.physical_sizes
@@ -153,11 +153,11 @@ def test_fill_undefined_data(periodic):
                             )
     assert topography.has_undefined_data
 
-    filled_topography = topography.fill_undefined_data(fill_value=-np.infty)
+    filled_topography = topography.fill_undefined_data(fill_value=-np.inf)
     mask = np.ma.getmask(topography.heights())
     nmask = np.logical_not(mask)
     assert (filled_topography[nmask] == topography[nmask]).all()
-    assert (filled_topography[mask] == - np.infty).all()
+    assert (filled_topography[mask] == - np.inf).all()
     assert not filled_topography.has_undefined_data
 
     assert filled_topography.physical_sizes == topography.physical_sizes
diff --git a/test/test_unit_conversion.py b/test/test_unit_conversion.py
index 098f4bc6b..c66b7afd5 100644
--- a/test/test_unit_conversion.py
+++ b/test/test_unit_conversion.py
@@ -25,10 +25,10 @@
 
 import numpy as np
 import pytest
-
 from NuMPI import MPI
 
-from SurfaceTopography.Support.UnitConversion import suggest_length_unit, suggest_length_unit_for_data
+from SurfaceTopography.Support.UnitConversion import (
+    suggest_length_unit, suggest_length_unit_for_data)
 
 pytestmark = pytest.mark.skipif(
     MPI.COMM_WORLD.Get_size() > 1,
@@ -56,8 +56,8 @@ def test_unit_outside_range():
 
 def test_nan_and_inf():
     assert suggest_length_unit_for_data('log', [1e-9, 10], 'm') == 'µm'
-    assert suggest_length_unit_for_data('log', [1e-9, 10, np.NaN], 'm') == 'µm'
-    assert suggest_length_unit_for_data('log', [1e-9, 10, np.Inf], 'm') == 'm'
-    assert suggest_length_unit_for_data('log', [1e-9, 10, np.Inf, np.Inf], 'm') == 'm'
+    assert suggest_length_unit_for_data('log', [1e-9, 10, np.nan], 'm') == 'µm'
+    assert suggest_length_unit_for_data('log', [1e-9, 10, np.inf], 'm') == 'm'
+    assert suggest_length_unit_for_data('log', [1e-9, 10, np.inf, np.inf], 'm') == 'm'
     with pytest.raises(ValueError):
-        suggest_length_unit_for_data('log', [np.NaN, np.NaN], 'm')
+        suggest_length_unit_for_data('log', [np.nan, np.nan], 'm')