diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index e116c483..f8dc74cd 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -16,16 +16,14 @@ jobs:
   source_check:
     name: source check
     runs-on: ubuntu-latest
-    strategy:
-      fail-fast: false
 
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
         with:
           fetch-depth: '0'
 
       - name: Set up Python 3.9
-        uses: actions/setup-python@v2
+        uses: actions/setup-python@v5
         with:
           python-version: 3.9
 
@@ -38,6 +36,10 @@ jobs:
         run: |
           python -m black --check --diff --color .
 
+      - name: black preview
+        run: |
+          python -m black --preview --diff --color .
+
       - name: isort check
         run: |
           python -m isort --check --diff --color .
@@ -58,60 +60,68 @@ jobs:
       matrix:
         cfg:
         - { os: ubuntu-latest, arch: x86_64 }
-        - { os: ubuntu-latest, arch: i686 }
         - { os: windows-latest, arch: AMD64 }
-        - { os: windows-latest, arch: x86 }
         - { os: macos-latest, arch: x86_64 }
         - { os: macos-latest, arch: arm64 }
         - { os: macos-latest, arch: universal2 }
 
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
         with:
           fetch-depth: '0'
 
       - name: Build wheels
-        uses: pypa/cibuildwheel@v2.16.2
+        uses: pypa/cibuildwheel@v2.17.0
         env:
           CIBW_ARCHS: ${{ matrix.cfg.arch }}
         with:
           output-dir: dist
 
-      - uses: actions/upload-artifact@v2
+      - uses: actions/upload-artifact@v3
         with:
           path: ./dist/*.whl
 
   build_sdist:
-    name: sdist on ${{ matrix.os }} with py ${{ matrix.python-version }}
+    name: sdist on ${{ matrix.os }} with py ${{ matrix.ver.py }} numpy${{ matrix.ver.np }} scipy${{ matrix.ver.sp }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, windows-latest, macos-latest]
-        python-version: ['3.8', '3.9', '3.10', '3.11', '3.12']
-
+        # https://github.com/scipy/oldest-supported-numpy/blob/main/setup.cfg
+        ver:
+          - {py: '3.8', np: '==1.20.0', sp: '==1.5.4'}
+          - {py: '3.9', np: '==1.20.0', sp: '==1.5.4'}
+          - {py: '3.10', np: '==1.21.6', sp: '==1.7.2'}
+          - {py: '3.11', np: '==1.23.2', sp: '==1.9.2'}
+          - {py: '3.12', np: '==1.26.2', sp: '==1.11.2'}
+          - {py: '3.12', np: '>=2.0.0rc1', sp: '>=1.13.0'}
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
         with:
           fetch-depth: '0'
 
-      - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v2
+      - name: Set up Python ${{ matrix.ver.py }}
+        uses: actions/setup-python@v5
         with:
-          python-version: ${{ matrix.python-version }}
+          python-version: ${{ matrix.ver.py }}
 
       - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install build "coveralls>=3.0.0"
+
+      - name: Install GSTools
         env:
           GSTOOLS_BUILD_PARALLEL: 1
         run: |
-          python -m pip install --upgrade pip
-          pip install build coveralls>=3.0.0
           pip install -v --editable .[test]
 
       - name: Run tests
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
         run: |
+          pip install "numpy${{ matrix.ver.np }}" "scipy${{ matrix.ver.sp }}"
           python -m pytest --cov gstools --cov-report term-missing -v tests/
           python -m coveralls --service=github
 
@@ -120,8 +130,8 @@ jobs:
           # PEP 517 package builder from pypa
           python -m build --sdist --outdir dist .
 
-      - uses: actions/upload-artifact@v2
-        if: matrix.os == 'ubuntu-latest' && matrix.python-version == '3.9'
+      - uses: actions/upload-artifact@v3
+        if: matrix.os == 'ubuntu-latest' && matrix.ver.py == '3.9'
         with:
           path: dist/*.tar.gz
 
@@ -130,7 +140,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-      - uses: actions/download-artifact@v2
+      - uses: actions/download-artifact@v3
         with:
           name: artifact
           path: dist
diff --git a/.gitignore b/.gitignore
index 872d3f40..bcdc980b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -117,6 +117,7 @@ src/gstools/_version.py
 # generated docs
 docs/source/examples/
 docs/source/api/
+docs/source/sg_execution_times.rst
 
 # other settings
 .vscode/
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 00513d23..61601c72 100755
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,25 @@
 
 All notable changes to **GSTools** will be documented in this file.
 
+## [Unreleased] - ?
+
+### Enhancements
+
+- added global variable `config.NUM_THREADS` to select number of threads for parallel computation ([#336](https://github.com/GeoStat-Framework/GSTools/pull/336))
+- prepare numpy 2 support ([#340](https://github.com/GeoStat-Framework/GSTools/pull/340))
+  - numpy 2.0.0rc1 for building extensions (for Python 3.9 and above)
+  - check multiple numpy and scipy versions in CI
+  - fixed minimal versions for numpy
+  - use `np.asarray` everywhere with `np.atleast_(n)d`
+  - fix long/longlong integer issue in cython on windows by always using 64bit integers
+
+### Bugfixes
+- build docs with latest sphinx version ([#340](https://github.com/GeoStat-Framework/GSTools/pull/340))
+
+### Changes
+- require pyvista 0.40 at least ([#340](https://github.com/GeoStat-Framework/GSTools/pull/340))
+
+
 ## [1.5.1] - Nifty Neon - 2023-11
 
 ### Enhancements
diff --git a/README.md b/README.md
index 69f1f3ec..6cb69901 100644
--- a/README.md
+++ b/README.md
@@ -345,7 +345,7 @@ in memory for immediate 3D plotting in Python.
 
 ## Requirements:
 
-- [NumPy >= 1.14.5](https://www.numpy.org)
+- [NumPy >= 1.20.0](https://www.numpy.org)
 - [SciPy >= 1.1.0](https://www.scipy.org/scipylib)
 - [hankel >= 1.0.0](https://github.com/steven-murray/hankel)
 - [emcee >= 3.0.0](https://github.com/dfm/emcee)
@@ -366,7 +366,7 @@ You can contact us via <info@geostat-framework.org>.
 
 ## License
 
-[LGPLv3][license_link] © 2018-2021
+[LGPLv3][license_link] © 2018-2024
 
 [pip_link]: https://pypi.org/project/gstools
 [conda_link]: https://docs.conda.io/en/latest/miniconda.html
diff --git a/docs/source/conf.py b/docs/source/conf.py
index eb5e5fae..e89928fc 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -259,10 +259,9 @@ def setup(app):
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
     "Python": ("https://docs.python.org/", None),
-    "NumPy": ("http://docs.scipy.org/doc/numpy/", None),
-    "SciPy": ("http://docs.scipy.org/doc/scipy/reference", None),
-    "matplotlib": ("http://matplotlib.org", None),
-    "Sphinx": ("http://www.sphinx-doc.org/en/stable/", None),
+    "NumPy": ("https://numpy.org/doc/stable/", None),
+    "SciPy": ("https://docs.scipy.org/doc/scipy/", None),
+    "matplotlib": ("https://matplotlib.org/stable/", None),
     "hankel": ("https://hankel.readthedocs.io/en/latest/", None),
     "emcee": ("https://emcee.readthedocs.io/en/latest/", None),
 }
diff --git a/docs/source/index.rst b/docs/source/index.rst
index d71494ac..ecad0583 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -414,7 +414,7 @@ in memory for immediate 3D plotting in Python.
 Requirements
 ============
 
-- `Numpy >= 1.14.5 <http://www.numpy.org>`_
+- `NumPy >= 1.20.0 <http://www.numpy.org>`_
 - `SciPy >= 1.1.0 <http://www.scipy.org>`_
 - `hankel >= 1.0.0 <https://github.com/steven-murray/hankel>`_
 - `emcee >= 3.0.0 <https://github.com/dfm/emcee>`_
diff --git a/examples/00_misc/00_tpl_stable.py b/examples/00_misc/00_tpl_stable.py
index 4d1a89d8..474b0f55 100644
--- a/examples/00_misc/00_tpl_stable.py
+++ b/examples/00_misc/00_tpl_stable.py
@@ -40,6 +40,7 @@
    \sigma^2_{\ell_{\mathrm{up}}} &=
    C\cdot\frac{\ell_{\mathrm{up}}^{2H}}{2H}
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/00_misc/01_export.py b/examples/00_misc/01_export.py
index ad754abf..e38294fe 100644
--- a/examples/00_misc/01_export.py
+++ b/examples/00_misc/01_export.py
@@ -7,6 +7,7 @@
 
 These can be viewed for example with `Paraview <https://www.paraview.org/>`__.
 """
+
 # sphinx_gallery_thumbnail_path = 'pics/paraview.png'
 import gstools as gs
 
diff --git a/examples/00_misc/02_check_rand_meth_sampling.py b/examples/00_misc/02_check_rand_meth_sampling.py
index 0b23a929..58d998b4 100644
--- a/examples/00_misc/02_check_rand_meth_sampling.py
+++ b/examples/00_misc/02_check_rand_meth_sampling.py
@@ -2,6 +2,7 @@
 Check Random Sampling
 ---------------------
 """
+
 import numpy as np
 from matplotlib import pyplot as plt
 from mpl_toolkits.mplot3d import Axes3D
diff --git a/examples/00_misc/04_herten.py b/examples/00_misc/04_herten.py
index d9e1b16e..1e1b8a23 100644
--- a/examples/00_misc/04_herten.py
+++ b/examples/00_misc/04_herten.py
@@ -30,6 +30,7 @@
 functions, since the only produce the ``herten_transmissivity.gz``
 and ``grid_dim_origin_spacing.txt``, which are already present.
 """
+
 import os
 
 import matplotlib.pyplot as plt
diff --git a/examples/00_misc/05_standalone_field.py b/examples/00_misc/05_standalone_field.py
index fb0c2082..e467f043 100644
--- a/examples/00_misc/05_standalone_field.py
+++ b/examples/00_misc/05_standalone_field.py
@@ -7,6 +7,7 @@
 In the following example we will produce 10000 random points in 4D with
 random values and plot them.
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/01_random_field/03_unstr_srf_export.py b/examples/01_random_field/03_unstr_srf_export.py
index e86e6f3e..94d00952 100644
--- a/examples/01_random_field/03_unstr_srf_export.py
+++ b/examples/01_random_field/03_unstr_srf_export.py
@@ -6,6 +6,7 @@
 Normally, such a grid would be read in, but we can simply generate one and
 then create a random field at those coordinates.
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/01_random_field/04_srf_merge.py b/examples/01_random_field/04_srf_merge.py
index 00aa6511..241ed079 100644
--- a/examples/01_random_field/04_srf_merge.py
+++ b/examples/01_random_field/04_srf_merge.py
@@ -6,6 +6,7 @@
 to merge two unstructured rectangular fields.
 
 """
+
 # sphinx_gallery_thumbnail_number = 2
 import numpy as np
 
diff --git a/examples/01_random_field/06_pyvista_support.py b/examples/01_random_field/06_pyvista_support.py
index 4cc371dd..29de8dd7 100644
--- a/examples/01_random_field/06_pyvista_support.py
+++ b/examples/01_random_field/06_pyvista_support.py
@@ -12,6 +12,7 @@
 The :any:`Field.mesh` method enables easy field creation on PyVista meshes
 used by the :any:`SRF` or :any:`Krige` class.
 """
+
 # sphinx_gallery_thumbnail_path = 'pics/GS_pyvista_cut.png'
 import pyvista as pv
 
@@ -21,8 +22,8 @@
 # We create a structured grid with PyVista containing 50 segments on all three
 # axes each with a length of 2 (whatever unit).
 
-dim, spacing = (50, 50, 50), (2, 2, 2)
-grid = pv.UniformGrid(dim, spacing)
+dims, spacing = (50, 50, 50), (2, 2, 2)
+grid = pv.ImageData(dimensions=dims, spacing=spacing)
 
 ###############################################################################
 # Now we set up the SRF class as always. We'll use an anisotropic model.
diff --git a/examples/02_cov_model/01_basic_methods.py b/examples/02_cov_model/01_basic_methods.py
index fb9bfe7a..4c97fba8 100755
--- a/examples/02_cov_model/01_basic_methods.py
+++ b/examples/02_cov_model/01_basic_methods.py
@@ -36,6 +36,7 @@
 correlation function as demonstrated in the introductory example.
 If one of the above functions is given, the others will be determined:
 """
+
 import gstools as gs
 
 model = gs.Exponential(dim=3, var=2.0, len_scale=10, nugget=0.5)
diff --git a/examples/02_cov_model/02_aniso_rotation.py b/examples/02_cov_model/02_aniso_rotation.py
index 2a8bac78..b7459e39 100755
--- a/examples/02_cov_model/02_aniso_rotation.py
+++ b/examples/02_cov_model/02_aniso_rotation.py
@@ -6,6 +6,7 @@
 represents the isotropic case for the model.
 Nevertheless, you can provide anisotropy ratios by:
 """
+
 import gstools as gs
 
 model = gs.Gaussian(dim=3, var=2.0, len_scale=10, anis=0.5)
diff --git a/examples/02_cov_model/03_spectral_methods.py b/examples/02_cov_model/03_spectral_methods.py
index 677811a9..61c7e49b 100755
--- a/examples/02_cov_model/03_spectral_methods.py
+++ b/examples/02_cov_model/03_spectral_methods.py
@@ -23,6 +23,7 @@
 
 You can access these methods by:
 """
+
 import gstools as gs
 
 model = gs.Gaussian(dim=3, var=2.0, len_scale=10)
diff --git a/examples/02_cov_model/04_different_scales.py b/examples/02_cov_model/04_different_scales.py
index cd6f4dee..0e2e1991 100755
--- a/examples/02_cov_model/04_different_scales.py
+++ b/examples/02_cov_model/04_different_scales.py
@@ -16,6 +16,7 @@
 
 You can access it by:
 """
+
 import gstools as gs
 
 model = gs.Stable(dim=3, var=2.0, len_scale=10)
diff --git a/examples/02_cov_model/05_additional_para.py b/examples/02_cov_model/05_additional_para.py
index 02507d96..3264cec4 100755
--- a/examples/02_cov_model/05_additional_para.py
+++ b/examples/02_cov_model/05_additional_para.py
@@ -10,6 +10,7 @@
 
 This leads to the so called **stable** covariance model and we can define it by
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/02_cov_model/06_fitting_para_ranges.py b/examples/02_cov_model/06_fitting_para_ranges.py
index bd3d1196..7c8f083e 100755
--- a/examples/02_cov_model/06_fitting_para_ranges.py
+++ b/examples/02_cov_model/06_fitting_para_ranges.py
@@ -6,6 +6,7 @@
 variogram data. In the following we will use the self defined stable model
 from a previous example.
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/03_variogram/00_fit_variogram.py b/examples/03_variogram/00_fit_variogram.py
index 9dc1b305..7334ed2c 100644
--- a/examples/03_variogram/00_fit_variogram.py
+++ b/examples/03_variogram/00_fit_variogram.py
@@ -2,6 +2,7 @@
 Fit Variogram
 -------------
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/03_variogram/01_find_best_model.py b/examples/03_variogram/01_find_best_model.py
index c0fc383f..eab031cf 100755
--- a/examples/03_variogram/01_find_best_model.py
+++ b/examples/03_variogram/01_find_best_model.py
@@ -2,6 +2,7 @@
 Finding the best fitting variogram model
 ----------------------------------------
 """
+
 import numpy as np
 from matplotlib import pyplot as plt
 
diff --git a/examples/03_variogram/02_multi_vario.py b/examples/03_variogram/02_multi_vario.py
index f893e452..71048849 100755
--- a/examples/03_variogram/02_multi_vario.py
+++ b/examples/03_variogram/02_multi_vario.py
@@ -5,6 +5,7 @@
 In this example, we demonstrate how to estimate a variogram from multiple
 fields on the same point-set that should have the same statistical properties.
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 
diff --git a/examples/03_variogram/03_directional_2d.py b/examples/03_variogram/03_directional_2d.py
index e5744cb2..460c7513 100755
--- a/examples/03_variogram/03_directional_2d.py
+++ b/examples/03_variogram/03_directional_2d.py
@@ -7,6 +7,7 @@
 
 Afterwards we will fit a model to this estimated variogram and show the result.
 """
+
 import numpy as np
 from matplotlib import pyplot as plt
 
diff --git a/examples/03_variogram/04_directional_3d.py b/examples/03_variogram/04_directional_3d.py
index 22ee7e9f..6a8b6ddf 100755
--- a/examples/03_variogram/04_directional_3d.py
+++ b/examples/03_variogram/04_directional_3d.py
@@ -7,6 +7,7 @@
 
 Afterwards we will fit a model to this estimated variogram and show the result.
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 from mpl_toolkits.mplot3d import Axes3D
diff --git a/examples/03_variogram/05_auto_fit_variogram.py b/examples/03_variogram/05_auto_fit_variogram.py
index 9c1cf2f9..2fcc7fbd 100644
--- a/examples/03_variogram/05_auto_fit_variogram.py
+++ b/examples/03_variogram/05_auto_fit_variogram.py
@@ -2,6 +2,7 @@
 Fit Variogram with automatic binning
 ------------------------------------
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/03_variogram/06_auto_bin_latlon.py b/examples/03_variogram/06_auto_bin_latlon.py
index 20f6f9bc..cc248ea1 100644
--- a/examples/03_variogram/06_auto_bin_latlon.py
+++ b/examples/03_variogram/06_auto_bin_latlon.py
@@ -8,6 +8,7 @@
 
 We use a data set from 20 meteo-stations choosen randomly.
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/04_vector_field/00_2d_vector_field.py b/examples/04_vector_field/00_2d_vector_field.py
index 6a6c71ac..2e722764 100644
--- a/examples/04_vector_field/00_2d_vector_field.py
+++ b/examples/04_vector_field/00_2d_vector_field.py
@@ -5,6 +5,7 @@
 As a first example we are going to generate a 2d vector field with a Gaussian
 covariance model on a structured grid:
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/04_vector_field/01_3d_vector_field.py b/examples/04_vector_field/01_3d_vector_field.py
index 13d76e87..5b1872bd 100755
--- a/examples/04_vector_field/01_3d_vector_field.py
+++ b/examples/04_vector_field/01_3d_vector_field.py
@@ -6,6 +6,7 @@
 Gaussian covariance model. The mesh on which we generate the field will be
 externally defined and it will be generated by PyVista.
 """
+
 # sphinx_gallery_thumbnail_path = 'pics/GS_3d_vector_field.png'
 import pyvista as pv
 
@@ -17,7 +18,7 @@
 ###############################################################################
 # create a uniform grid with PyVista
 dims, spacing, origin = (40, 30, 10), (1, 1, 1), (-10, 0, 0)
-mesh = pv.UniformGrid(dims=dims, spacing=spacing, origin=origin)
+mesh = pv.ImageData(dimensions=dims, spacing=spacing, origin=origin)
 
 ###############################################################################
 # create an incompressible random 3d velocity field on the given mesh
diff --git a/examples/05_kriging/00_simple_kriging.py b/examples/05_kriging/00_simple_kriging.py
index bb83a32c..1a245b6c 100755
--- a/examples/05_kriging/00_simple_kriging.py
+++ b/examples/05_kriging/00_simple_kriging.py
@@ -27,6 +27,7 @@
 The mean of the field has to be given beforehand.
 
 """
+
 import numpy as np
 
 from gstools import Gaussian, krige
diff --git a/examples/05_kriging/01_ordinary_kriging.py b/examples/05_kriging/01_ordinary_kriging.py
index 3bc54bd9..d26254ef 100644
--- a/examples/05_kriging/01_ordinary_kriging.py
+++ b/examples/05_kriging/01_ordinary_kriging.py
@@ -27,6 +27,7 @@
 Here we use ordinary kriging in 1D (for plotting reasons) with 5 given observations/conditions.
 The estimated mean can be accessed by ``krig.mean``.
 """
+
 import numpy as np
 
 from gstools import Gaussian, krige
diff --git a/examples/05_kriging/02_pykrige_interface.py b/examples/05_kriging/02_pykrige_interface.py
index bb166d29..a6fbf03e 100755
--- a/examples/05_kriging/02_pykrige_interface.py
+++ b/examples/05_kriging/02_pykrige_interface.py
@@ -11,6 +11,7 @@
 To demonstrate the general workflow, we compare ordinary kriging of PyKrige
 with the corresponding GSTools routine in 2D:
 """
+
 import numpy as np
 from matplotlib import pyplot as plt
 from pykrige.ok import OrdinaryKriging
diff --git a/examples/05_kriging/03_compare_kriging.py b/examples/05_kriging/03_compare_kriging.py
index 12f5c43b..463faa0a 100755
--- a/examples/05_kriging/03_compare_kriging.py
+++ b/examples/05_kriging/03_compare_kriging.py
@@ -2,6 +2,7 @@
 Compare Kriging
 ---------------
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 
diff --git a/examples/05_kriging/04_extdrift_kriging.py b/examples/05_kriging/04_extdrift_kriging.py
index b9bbb54a..2e6a168a 100755
--- a/examples/05_kriging/04_extdrift_kriging.py
+++ b/examples/05_kriging/04_extdrift_kriging.py
@@ -2,6 +2,7 @@
 External Drift Kriging
 ----------------------
 """
+
 import numpy as np
 
 from gstools import SRF, Gaussian, krige
diff --git a/examples/05_kriging/05_universal_kriging.py b/examples/05_kriging/05_universal_kriging.py
index d3cf2880..5501694a 100755
--- a/examples/05_kriging/05_universal_kriging.py
+++ b/examples/05_kriging/05_universal_kriging.py
@@ -13,6 +13,7 @@
 To access only the estimated mean/drift, we provide a switch `only_mean`
 in the call routine.
 """
+
 import numpy as np
 
 from gstools import SRF, Gaussian, krige
diff --git a/examples/05_kriging/06_detrended_kriging.py b/examples/05_kriging/06_detrended_kriging.py
index b68e6502..6d20cf1d 100755
--- a/examples/05_kriging/06_detrended_kriging.py
+++ b/examples/05_kriging/06_detrended_kriging.py
@@ -2,6 +2,7 @@
 Detrended Kriging
 -----------------
 """
+
 import numpy as np
 
 from gstools import SRF, Gaussian, krige
diff --git a/examples/05_kriging/07_detrended_ordinary_kriging.py b/examples/05_kriging/07_detrended_ordinary_kriging.py
index 9a59fc33..81d01744 100755
--- a/examples/05_kriging/07_detrended_ordinary_kriging.py
+++ b/examples/05_kriging/07_detrended_ordinary_kriging.py
@@ -2,6 +2,7 @@
 Detrended Ordinary Kriging
 --------------------------
 """
+
 import numpy as np
 
 from gstools import SRF, Gaussian, krige
diff --git a/examples/05_kriging/09_pseudo_inverse.py b/examples/05_kriging/09_pseudo_inverse.py
index 9920a956..7615d888 100755
--- a/examples/05_kriging/09_pseudo_inverse.py
+++ b/examples/05_kriging/09_pseudo_inverse.py
@@ -17,6 +17,7 @@
 In the following we have two different values at the same location.
 The resulting kriging field will hold the average at this point.
 """
+
 import numpy as np
 
 from gstools import Gaussian, krige
diff --git a/examples/06_conditioned_fields/00_condition_ensemble.py b/examples/06_conditioned_fields/00_condition_ensemble.py
index 716edd20..5cc07eed 100644
--- a/examples/06_conditioned_fields/00_condition_ensemble.py
+++ b/examples/06_conditioned_fields/00_condition_ensemble.py
@@ -6,6 +6,7 @@
 with 5 given observations/conditions,
 to generate an ensemble of conditioned random fields.
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 
diff --git a/examples/06_conditioned_fields/01_2D_condition_ensemble.py b/examples/06_conditioned_fields/01_2D_condition_ensemble.py
index 4a03e66b..81f51464 100644
--- a/examples/06_conditioned_fields/01_2D_condition_ensemble.py
+++ b/examples/06_conditioned_fields/01_2D_condition_ensemble.py
@@ -4,6 +4,7 @@
 
 Let's create an ensemble of conditioned random fields in 2D.
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 
diff --git a/examples/07_transformations/00_log_normal.py b/examples/07_transformations/00_log_normal.py
index 3d2f4532..d44c1627 100755
--- a/examples/07_transformations/00_log_normal.py
+++ b/examples/07_transformations/00_log_normal.py
@@ -6,6 +6,7 @@
 
 See :any:`transform.normal_to_lognormal`
 """
+
 import gstools as gs
 
 # structured field with a size of 100x100 and a grid-size of 1x1
diff --git a/examples/07_transformations/01_binary.py b/examples/07_transformations/01_binary.py
index a403abb7..125e29d0 100755
--- a/examples/07_transformations/01_binary.py
+++ b/examples/07_transformations/01_binary.py
@@ -8,6 +8,7 @@
 
 See :any:`transform.binary`
 """
+
 import gstools as gs
 
 # structured field with a size of 100x100 and a grid-size of 1x1
diff --git a/examples/07_transformations/02_discrete.py b/examples/07_transformations/02_discrete.py
index 268cc960..48f67a2d 100755
--- a/examples/07_transformations/02_discrete.py
+++ b/examples/07_transformations/02_discrete.py
@@ -9,6 +9,7 @@
 
 See :any:`transform.discrete`
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/07_transformations/03_zinn_harvey.py b/examples/07_transformations/03_zinn_harvey.py
index 01e65790..fad1fb64 100755
--- a/examples/07_transformations/03_zinn_harvey.py
+++ b/examples/07_transformations/03_zinn_harvey.py
@@ -9,6 +9,7 @@
 
 See :any:`transform.zinnharvey`
 """
+
 import gstools as gs
 
 # structured field with a size of 100x100 and a grid-size of 1x1
diff --git a/examples/07_transformations/04_bimodal.py b/examples/07_transformations/04_bimodal.py
index 8234a486..4dd6fb29 100755
--- a/examples/07_transformations/04_bimodal.py
+++ b/examples/07_transformations/04_bimodal.py
@@ -11,6 +11,7 @@
 
 See: :any:`transform.normal_to_arcsin` and :any:`transform.normal_to_uquad`
 """
+
 import gstools as gs
 
 # structured field with a size of 100x100 and a grid-size of 1x1
diff --git a/examples/07_transformations/05_combinations.py b/examples/07_transformations/05_combinations.py
index fb27c5e6..1fbe367e 100755
--- a/examples/07_transformations/05_combinations.py
+++ b/examples/07_transformations/05_combinations.py
@@ -15,6 +15,7 @@
 
 If you don't specify `field` and `store` everything happens inplace.
 """
+
 # sphinx_gallery_thumbnail_number = 1
 import gstools as gs
 
diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py
index b7b3748a..5b1a6fca 100755
--- a/examples/08_geo_coordinates/00_field_generation.py
+++ b/examples/08_geo_coordinates/00_field_generation.py
@@ -16,6 +16,7 @@
 To generate the field, we simply pass ``(lat, lon)`` as the position tuple
 to the :any:`SRF` class.
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py
index 44b9f28c..3c17fb7e 100755
--- a/examples/08_geo_coordinates/01_dwd_krige.py
+++ b/examples/08_geo_coordinates/01_dwd_krige.py
@@ -15,6 +15,7 @@
 In order to keep the number of dependecies low, the calls of both functions
 shown beneath are commented out.
 """
+
 # sphinx_gallery_thumbnail_number = 2
 import matplotlib.pyplot as plt
 import numpy as np
diff --git a/examples/09_spatio_temporal/03_geographic_coordinates.py b/examples/09_spatio_temporal/03_geographic_coordinates.py
index 56684a94..b1cfbff6 100644
--- a/examples/09_spatio_temporal/03_geographic_coordinates.py
+++ b/examples/09_spatio_temporal/03_geographic_coordinates.py
@@ -18,6 +18,7 @@
 
 We will set a spatial length-scale of `1000` and a time length-scale of `100` days.
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/10_normalizer/00_lognormal_kriging.py b/examples/10_normalizer/00_lognormal_kriging.py
index b7bf4a8a..9880bc37 100644
--- a/examples/10_normalizer/00_lognormal_kriging.py
+++ b/examples/10_normalizer/00_lognormal_kriging.py
@@ -14,6 +14,7 @@
 
 In this example we will use ordinary kriging.
 """
+
 import numpy as np
 
 import gstools as gs
diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py
index 80a49569..71ad1385 100644
--- a/examples/10_normalizer/01_auto_fit.py
+++ b/examples/10_normalizer/01_auto_fit.py
@@ -18,6 +18,7 @@
 We will generate the "original" field on a 60x60 mesh, from which we will take
 samples in order to pretend a situation of data-scarcity.
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 
diff --git a/examples/10_normalizer/02_compare.py b/examples/10_normalizer/02_compare.py
index 24f7c9f8..2dd74488 100644
--- a/examples/10_normalizer/02_compare.py
+++ b/examples/10_normalizer/02_compare.py
@@ -6,6 +6,7 @@
 
 But first, we define a convenience routine and make some imports as always.
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 
diff --git a/pyproject.toml b/pyproject.toml
index 70b7d0a6..17bcb246 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -2,8 +2,9 @@
 requires = [
     "setuptools>=64",
     "setuptools_scm>=7",
-    "oldest-supported-numpy",
-    "Cython>=3.0",
+    "numpy>=2.0.0rc1,<2.3; python_version >= '3.9'",
+    "oldest-supported-numpy; python_version < '3.9'",
+    "Cython>=3.0.10,<3.1.0",
     "extension-helpers>=1",
 ]
 build-backend = "setuptools.build_meta"
@@ -48,7 +49,7 @@ dependencies = [
     "emcee>=3.0.0",
     "hankel>=1.0.0",
     "meshio>=5.1.0",
-    "numpy>=1.14.5",
+    "numpy>=1.20.0",
     "pyevtk>=1.1.1",
     "scipy>=1.1.0",
 ]
@@ -60,20 +61,20 @@ doc = [
     "meshzoo>=0.7",
     "numpydoc>=1.1",
     "pykrige>=1.5,<2",
-    "pyvista>=0.29",
-    "sphinx>=4",
+    "pyvista>=0.40",
+    "sphinx>=7",
     "sphinx-gallery>=0.8",
-    "sphinx-rtd-theme>=1,<1.1",
+    "sphinx-rtd-theme>=2",
     "sphinxcontrib-youtube>=1.1",
 ]
 plotting = [
     "matplotlib>=3",
-    "pyvista>=0.29",
+    "pyvista>=0.40",
 ]
 rust = ["gstools_core>=0.2.0,<1"]
 test = ["pytest-cov>=3"]
 lint = [
-    "black",
+    "black>=24",
     "pylint",
     "isort[colors]",
     "cython-lint",
@@ -161,11 +162,9 @@ target-version = [
 # Switch to using build
 build-frontend = "build"
 # Disable building PyPy wheels on all platforms, 32bit for py3.10/11/12, musllinux builds, py3.6/7
-skip = ["cp36-*", "cp37-*", "pp*", "cp31*-win32", "cp31*-manylinux_i686", "*-musllinux_*"]
+skip = ["cp36-*", "cp37-*", "pp*", "*-win32", "*-manylinux_i686", "*-musllinux_*"]
 # Run the package tests using `pytest`
 test-extras = "test"
 test-command = "pytest -v {package}/tests"
 # Skip trying to test arm64 builds on Intel Macs
 test-skip = "*-macosx_arm64 *-macosx_universal2:arm64"
-# no wheels for linux-32bit anymore for numpy>=1.22
-environment = "PIP_PREFER_BINARY=1"
diff --git a/setup.py b/setup.py
index f465ffd4..b27548a9 100644
--- a/setup.py
+++ b/setup.py
@@ -1,4 +1,5 @@
 """GSTools: A geostatistical toolbox."""
+
 import os
 
 import numpy as np
diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py
index 6d62d558..11e63a2b 100644
--- a/src/gstools/__init__.py
+++ b/src/gstools/__init__.py
@@ -129,6 +129,7 @@
    DEGREE_SCALE
    RADIAN_SCALE
 """
+
 # Hooray!
 from gstools import (
     config,
diff --git a/src/gstools/config.py b/src/gstools/config.py
index f0a87bb6..24ce20c7 100644
--- a/src/gstools/config.py
+++ b/src/gstools/config.py
@@ -4,6 +4,7 @@
 .. currentmodule:: gstools.config
 
 """
+
 NUM_THREADS = None
 
 # pylint: disable=W0611
diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py
index be530075..23e19881 100644
--- a/src/gstools/covmodel/base.py
+++ b/src/gstools/covmodel/base.py
@@ -8,6 +8,7 @@
 .. autosummary::
    CovModel
 """
+
 # pylint: disable=C0103, R0201, E1101, C0302, W0613
 import copy
 
diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py
index dc2d5a3a..b420be70 100755
--- a/src/gstools/covmodel/fit.py
+++ b/src/gstools/covmodel/fit.py
@@ -8,6 +8,7 @@
 .. autosummary::
    fit_variogram
 """
+
 # pylint: disable=C0103, W0632
 import numpy as np
 from scipy.optimize import curve_fit
@@ -308,9 +309,11 @@ def _pre_init_guess(model, init_guess, mean_x=1.0, mean_y=1.0):
     for opt in model.opt_arg:
         init_guess.setdefault(
             opt,
-            default_arg_from_bounds(bnd[opt])
-            if default
-            else getattr(model, opt),
+            (
+                default_arg_from_bounds(bnd[opt])
+                if default
+                else getattr(model, opt)
+            ),
         )
     # convert all init guesses to float (except "anis")
     for arg in model.iso_arg:
diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py
index a2aaf6a3..ab94e279 100644
--- a/src/gstools/covmodel/models.py
+++ b/src/gstools/covmodel/models.py
@@ -20,6 +20,7 @@
    SuperSpherical
    JBessel
 """
+
 # pylint: disable=C0103, E1101, R0201
 import warnings
 
diff --git a/src/gstools/covmodel/plot.py b/src/gstools/covmodel/plot.py
index 334ccbe2..32148c14 100644
--- a/src/gstools/covmodel/plot.py
+++ b/src/gstools/covmodel/plot.py
@@ -24,6 +24,7 @@
    plot_spectral_density
    plot_spectral_rad_pdf
 """
+
 # pylint: disable=C0103, C0415, E1130
 import numpy as np
 
diff --git a/src/gstools/covmodel/tpl_models.py b/src/gstools/covmodel/tpl_models.py
index 5e980760..b728e7b9 100644
--- a/src/gstools/covmodel/tpl_models.py
+++ b/src/gstools/covmodel/tpl_models.py
@@ -11,6 +11,7 @@
    TPLStable
    TPLSimple
 """
+
 # pylint: disable=C0103, E1101
 import warnings
 
diff --git a/src/gstools/field/base.py b/src/gstools/field/base.py
index bb514141..2006b858 100755
--- a/src/gstools/field/base.py
+++ b/src/gstools/field/base.py
@@ -8,6 +8,7 @@
 .. autosummary::
    Field
 """
+
 # pylint: disable=C0103, C0415
 from collections.abc import Iterable
 from copy import copy
@@ -362,7 +363,7 @@ def transform(
         ----------
         method : :class:`str`
             Method to use.
-            See :any:`gstools.transform` for available transformations.
+            See :py:mod:`gstools.transform` for available transformations.
         field : :class:`str`, optional
             Name of field to be transformed. The default is "field".
         store : :class:`str` or :class:`bool`, optional
diff --git a/src/gstools/field/cond_srf.py b/src/gstools/field/cond_srf.py
index 385216ba..c3e03fe2 100644
--- a/src/gstools/field/cond_srf.py
+++ b/src/gstools/field/cond_srf.py
@@ -8,6 +8,7 @@
 .. autosummary::
    CondSRF
 """
+
 # pylint: disable=C0103, W0231, W0221, W0222, E1102
 
 import numpy as np
diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py
index 571c9733..5beab10d 100644
--- a/src/gstools/field/generator.py
+++ b/src/gstools/field/generator.py
@@ -12,6 +12,7 @@
    RandMeth
    IncomprRandMeth
 """
+
 # pylint: disable=C0103, W0222, C0412, W0231
 import warnings
 from abc import ABC, abstractmethod
diff --git a/src/gstools/field/plot.py b/src/gstools/field/plot.py
index d06a22a1..ab28c974 100644
--- a/src/gstools/field/plot.py
+++ b/src/gstools/field/plot.py
@@ -9,6 +9,7 @@
    plot_field
    plot_vec_field
 """
+
 # pylint: disable=C0103, W0613, E1101
 import numpy as np
 from scipy import interpolate as inter
diff --git a/src/gstools/field/srf.py b/src/gstools/field/srf.py
index a8a1e575..d88e46c0 100644
--- a/src/gstools/field/srf.py
+++ b/src/gstools/field/srf.py
@@ -8,6 +8,7 @@
 .. autosummary::
    SRF
 """
+
 # pylint: disable=C0103, W0221, E1102
 
 import numpy as np
diff --git a/src/gstools/field/tools.py b/src/gstools/field/tools.py
index 4b128d8a..dfa2e3c6 100644
--- a/src/gstools/field/tools.py
+++ b/src/gstools/field/tools.py
@@ -10,6 +10,7 @@
    to_vtk_helper
    generate_on_mesh
 """
+
 # pylint: disable=W0212, C0415
 import meshio
 import numpy as np
diff --git a/src/gstools/field/upscaling.py b/src/gstools/field/upscaling.py
index 6dcbaa7b..857bfc45 100644
--- a/src/gstools/field/upscaling.py
+++ b/src/gstools/field/upscaling.py
@@ -11,6 +11,7 @@
    var_coarse_graining
    var_no_scaling
 """
+
 # pylint: disable=W0613
 import warnings
 
diff --git a/src/gstools/krige/__init__.py b/src/gstools/krige/__init__.py
index bb3ef699..66d03246 100644
--- a/src/gstools/krige/__init__.py
+++ b/src/gstools/krige/__init__.py
@@ -16,6 +16,7 @@
    ExtDrift
    Detrended
 """
+
 from gstools.krige.base import Krige
 from gstools.krige.methods import (
     Detrended,
diff --git a/src/gstools/krige/base.py b/src/gstools/krige/base.py
index 336f4cc0..49a4f62f 100755
--- a/src/gstools/krige/base.py
+++ b/src/gstools/krige/base.py
@@ -8,6 +8,7 @@
 .. autosummary::
    Krige
 """
+
 # pylint: disable=C0103, W0221, E1102, R0201, C0412
 import collections
 
@@ -383,9 +384,7 @@ def _pre_ext_drift(self, pnt_cnt, ext_drift=None, set_cond=False):
             the drift values at the given positions
         """
         if ext_drift is not None:
-            ext_drift = np.array(
-                ext_drift, dtype=np.double, ndmin=2, copy=False
-            )
+            ext_drift = np.atleast_2d(np.asarray(ext_drift, dtype=np.double))
             if ext_drift.size == 0:  # treat empty array as no ext_drift
                 return np.array([])
             if set_cond:
diff --git a/src/gstools/krige/methods.py b/src/gstools/krige/methods.py
index b258a02d..19ffed56 100644
--- a/src/gstools/krige/methods.py
+++ b/src/gstools/krige/methods.py
@@ -12,6 +12,7 @@
    ExtDrift
    Detrended
 """
+
 # pylint: disable=C0103
 from gstools.krige.base import Krige
 
diff --git a/src/gstools/krige/tools.py b/src/gstools/krige/tools.py
index e3112ae4..62926595 100644
--- a/src/gstools/krige/tools.py
+++ b/src/gstools/krige/tools.py
@@ -9,6 +9,7 @@
    set_condition
    get_drift_functions
 """
+
 # pylint: disable=C0103
 from itertools import combinations_with_replacement
 
diff --git a/src/gstools/normalizer/base.py b/src/gstools/normalizer/base.py
index 0072d6a9..4a8477c6 100644
--- a/src/gstools/normalizer/base.py
+++ b/src/gstools/normalizer/base.py
@@ -8,6 +8,7 @@
 .. autosummary::
    Normalizer
 """
+
 # pylint: disable=R0201
 import warnings
 
diff --git a/src/gstools/normalizer/methods.py b/src/gstools/normalizer/methods.py
index f66cfd19..a46dc230 100644
--- a/src/gstools/normalizer/methods.py
+++ b/src/gstools/normalizer/methods.py
@@ -13,6 +13,7 @@
    Modulus
    Manly
 """
+
 # pylint: disable=E1101
 import numpy as np
 
diff --git a/src/gstools/normalizer/tools.py b/src/gstools/normalizer/tools.py
index 7a5237b7..3e395d29 100644
--- a/src/gstools/normalizer/tools.py
+++ b/src/gstools/normalizer/tools.py
@@ -9,6 +9,7 @@
    apply_mean_norm_trend
    remove_trend_norm_mean
 """
+
 import numpy as np
 
 from gstools.normalizer.base import Normalizer
diff --git a/src/gstools/random/rng.py b/src/gstools/random/rng.py
index 4f698557..ad78a0c6 100644
--- a/src/gstools/random/rng.py
+++ b/src/gstools/random/rng.py
@@ -8,6 +8,7 @@
 .. autosummary::
    RNG
 """
+
 # pylint: disable=E1101
 import emcee as mc
 import numpy as np
diff --git a/src/gstools/tools/export.py b/src/gstools/tools/export.py
index 3e522c92..38254ceb 100644
--- a/src/gstools/tools/export.py
+++ b/src/gstools/tools/export.py
@@ -13,6 +13,7 @@
    to_vtk_structured
    to_vtk_unstructured
 """
+
 # pylint: disable=C0103, E1101
 import numpy as np
 from pyevtk.hl import gridToVTK, pointsToVTK
diff --git a/src/gstools/tools/geometric.py b/src/gstools/tools/geometric.py
index 7f2ea10e..55408965 100644
--- a/src/gstools/tools/geometric.py
+++ b/src/gstools/tools/geometric.py
@@ -29,6 +29,7 @@
    chordal_to_great_circle
    great_circle_to_chordal
 """
+
 # pylint: disable=C0103
 import numpy as np
 
@@ -383,7 +384,7 @@ def generate_st_grid(pos, time, mesh_type="unstructured"):
     if mesh_type != "unstructured":
         pos = generate_grid(pos)
     else:
-        pos = np.array(pos, dtype=np.double, ndmin=2, copy=False)
+        pos = np.atleast_2d(np.asarray(pos, dtype=np.double))
     out = [np.repeat(p.reshape(-1), np.size(time)) for p in pos]
     out.append(np.tile(time, np.size(pos[0])))
     return np.asarray(out, dtype=np.double)
@@ -552,7 +553,7 @@ def format_unstruct_pos_shape(pos, shape, check_stacked_shape=False):
     # now we try to be smart
     pre_len = len(np.atleast_1d(pos))
     # care about 1D: pos can be given as 1D array here -> convert to 2D array
-    pos = np.array(pos, dtype=np.double, ndmin=2, copy=False)
+    pos = np.atleast_2d(np.asarray(pos, dtype=np.double))
     post_len = len(pos)
     # first array dimension should be spatial dimension (1D is special case)
     dim = post_len if pre_len == post_len else 1
@@ -606,7 +607,7 @@ def ang2dir(angles, dtype=np.double, dim=None):
         the array of direction vectors
     """
     pre_dim = np.asanyarray(angles).ndim
-    angles = np.array(angles, ndmin=2, dtype=dtype, copy=False)
+    angles = np.atleast_2d(np.asarray(angles, dtype=dtype))
     if len(angles.shape) > 2:
         raise ValueError(f"Can't interpret angles array {angles}")
     dim = angles.shape[1] + 1 if dim is None else dim
diff --git a/src/gstools/tools/misc.py b/src/gstools/tools/misc.py
index 068200b7..aaba1501 100755
--- a/src/gstools/tools/misc.py
+++ b/src/gstools/tools/misc.py
@@ -10,6 +10,7 @@
    list_format
    eval_func
 """
+
 # pylint: disable=C0103, C0415
 import numpy as np
 
diff --git a/src/gstools/tools/special.py b/src/gstools/tools/special.py
index c48887ea..1457b736 100644
--- a/src/gstools/tools/special.py
+++ b/src/gstools/tools/special.py
@@ -14,6 +14,7 @@
    tpl_exp_spec_dens
    tpl_gau_spec_dens
 """
+
 # pylint: disable=C0103, E1101
 import numpy as np
 from scipy import special as sps
diff --git a/src/gstools/transform/array.py b/src/gstools/transform/array.py
index 30c92bb6..87564edf 100644
--- a/src/gstools/transform/array.py
+++ b/src/gstools/transform/array.py
@@ -18,6 +18,7 @@
    array_to_arcsin
    array_to_uquad
 """
+
 # pylint: disable=C0103, C0123, R0911
 from warnings import warn
 
diff --git a/src/gstools/transform/field.py b/src/gstools/transform/field.py
index 5d204407..4a281564 100644
--- a/src/gstools/transform/field.py
+++ b/src/gstools/transform/field.py
@@ -26,6 +26,7 @@
    normal_to_arcsin
    normal_to_uquad
 """
+
 # pylint: disable=C0103, C0123, R0911, R1735
 import numpy as np
 
@@ -111,7 +112,7 @@ def apply(fld, method, field="field", store=True, process=False, **kwargs):
         Field class containing a generated field.
     method : :class:`str`
         Method to use.
-        See :any:`gstools.transform` for available transformations.
+        See :py:mod:`gstools.transform` for available transformations.
     field : :class:`str`, optional
         Name of field to be transformed. The default is "field".
     store : :class:`str` or :class:`bool`, optional
diff --git a/src/gstools/variogram/binning.py b/src/gstools/variogram/binning.py
index e8e42f38..86d4fdc2 100644
--- a/src/gstools/variogram/binning.py
+++ b/src/gstools/variogram/binning.py
@@ -8,6 +8,7 @@
 .. autosummary::
    standard_bins
 """
+
 import numpy as np
 
 from gstools.tools import RADIAN_SCALE
diff --git a/src/gstools/variogram/estimator.pyx b/src/gstools/variogram/estimator.pyx
index ab418e32..e00824be 100644
--- a/src/gstools/variogram/estimator.pyx
+++ b/src/gstools/variogram/estimator.pyx
@@ -116,7 +116,7 @@ ctypedef double (*_estimator_func)(const double) nogil
 
 cdef inline void normalization_matheron(
     double[:] variogram,
-    long[:] counts,
+    np.int64_t[:] counts,
 ):
     cdef int i
     for i in range(variogram.shape[0]):
@@ -125,10 +125,10 @@ cdef inline void normalization_matheron(
 
 cdef inline void normalization_cressie(
     double[:] variogram,
-    long[:] counts,
+    np.int64_t[:] counts,
 ):
     cdef int i
-    cdef long cnt
+    cdef np.int64_t cnt
     for i in range(variogram.shape[0]):
         # avoid division by zero
         cnt = max(counts[i], 1)
@@ -139,12 +139,12 @@ cdef inline void normalization_cressie(
 
 ctypedef void (*_normalization_func)(
     double[:],
-    long[:],
+    np.int64_t[:],
 )
 
 cdef inline void normalization_matheron_vec(
     double[:, :] variogram,
-    long[:, :] counts,
+    np.int64_t[:, :] counts,
 ):
     cdef int d
     for d in range(variogram.shape[0]):
@@ -152,7 +152,7 @@ cdef inline void normalization_matheron_vec(
 
 cdef inline void normalization_cressie_vec(
     double[:, :] variogram,
-    long[:, :] counts,
+    np.int64_t[:, :] counts,
 ):
     cdef int d
     for d in range(variogram.shape[0]):
@@ -160,7 +160,7 @@ cdef inline void normalization_cressie_vec(
 
 ctypedef void (*_normalization_func_vec)(
     double[:, :],
-    long[:, :],
+    np.int64_t[:, :],
 )
 
 cdef _estimator_func choose_estimator_func(str estimator_type):
@@ -221,7 +221,7 @@ def directional(
     cdef int f_max = f.shape[0]
 
     cdef double[:, :] variogram = np.zeros((d_max, len(bin_edges)-1))
-    cdef long[:, :] counts = np.zeros((d_max, len(bin_edges)-1), dtype=long)
+    cdef np.int64_t[:, :] counts = np.zeros((d_max, len(bin_edges)-1), dtype=np.int64)
     cdef int i, j, k, m, d
     cdef double dist
 
@@ -287,7 +287,7 @@ def unstructured(
     cdef int f_max = f.shape[0]
 
     cdef double[:] variogram = np.zeros(len(bin_edges)-1)
-    cdef long[:] counts = np.zeros(len(bin_edges)-1, dtype=long)
+    cdef np.int64_t[:] counts = np.zeros(len(bin_edges)-1, dtype=np.int64)
     cdef int i, j, k, m
     cdef double dist
 
@@ -324,7 +324,7 @@ def structured(
     cdef int k_max = i_max + 1
 
     cdef double[:] variogram = np.zeros(k_max)
-    cdef long[:] counts = np.zeros(k_max, dtype=long)
+    cdef np.int64_t[:] counts = np.zeros(k_max, dtype=np.int64)
     cdef int i, j, k
 
     cdef int num_threads_c = set_num_threads(num_threads)
@@ -356,7 +356,7 @@ def ma_structured(
     cdef int k_max = i_max + 1
 
     cdef double[:] variogram = np.zeros(k_max)
-    cdef long[:] counts = np.zeros(k_max, dtype=long)
+    cdef np.int64_t[:] counts = np.zeros(k_max, dtype=np.int64)
     cdef int i, j, k
 
     cdef int num_threads_c = set_num_threads(num_threads)
diff --git a/src/gstools/variogram/variogram.py b/src/gstools/variogram/variogram.py
index 31e29191..afcf336f 100644
--- a/src/gstools/variogram/variogram.py
+++ b/src/gstools/variogram/variogram.py
@@ -9,6 +9,7 @@
    vario_estimate
    vario_estimate_axis
 """
+
 # pylint: disable=C0412
 import numpy as np
 
@@ -262,7 +263,7 @@ def vario_estimate(
            John Wiley & Sons. (2007)
     """
     if bin_edges is not None:
-        bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double, copy=False)
+        bin_edges = np.atleast_1d(np.asarray(bin_edges, dtype=np.double))
         bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0
     # allow multiple fields at same positions (ndmin=2: first axis -> field ID)
     # need to convert to ma.array, since list of ma.array is not recognised
@@ -315,7 +316,7 @@ def vario_estimate(
     # set directions
     dir_no = 0
     if direction is not None and dim > 1:
-        direction = np.array(direction, ndmin=2, dtype=np.double, copy=False)
+        direction = np.atleast_2d(np.asarray(direction, dtype=np.double))
         if len(direction.shape) > 2:
             raise ValueError(f"Can't interpret directions: {direction}")
         if direction.shape[1] != dim:
@@ -473,7 +474,7 @@ def vario_estimate_axis(
         if not config.USE_RUST:
             mask = np.asarray(mask, dtype=np.int32)
     else:
-        field = np.array(field, ndmin=1, dtype=np.double, copy=False)
+        field = np.atleast_1d(np.asarray(field, dtype=np.double))
         missing_mask = None  # free space
 
     axis_to_swap = AXIS_DIR[direction] if direction in AXIS else int(direction)
diff --git a/tests/test_condition.py b/tests/test_condition.py
index e6b0c0cc..8d5d0535 100644
--- a/tests/test_condition.py
+++ b/tests/test_condition.py
@@ -1,4 +1,5 @@
 """This is the unittest of CondSRF class."""
+
 import unittest
 from copy import copy
 
diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py
index b4344cc2..a2729dd6 100644
--- a/tests/test_covmodel.py
+++ b/tests/test_covmodel.py
@@ -1,6 +1,7 @@
 """
 This is the unittest of CovModel class.
 """
+
 import unittest
 
 import numpy as np
diff --git a/tests/test_export.py b/tests/test_export.py
index 9bd6dd41..b32898f4 100644
--- a/tests/test_export.py
+++ b/tests/test_export.py
@@ -1,5 +1,5 @@
-"""Test the PyVista/VTK export methods
-"""
+"""Test the PyVista/VTK export methods"""
+
 import os
 import shutil
 import tempfile
diff --git a/tests/test_transform.py b/tests/test_transform.py
index eb2980e2..abc5505b 100644
--- a/tests/test_transform.py
+++ b/tests/test_transform.py
@@ -1,4 +1,5 @@
 """This is the unittest of the transform submodule."""
+
 import unittest
 
 import numpy as np