diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index e2878e2f..7271ad70 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -17,7 +17,7 @@ jobs:
     strategy:
       matrix:
         os: ["ubuntu-latest"]
-        python-version: ["3.7", "3.8", "3.9"]
+        python-version: ["3.8", "3.9", "3.10"]
 
 
     # Run all shells using bash (including Windows)
diff --git a/.readthedocs.yml b/.readthedocs.yml
index 07df20ac..330ebc12 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yml
@@ -14,7 +14,7 @@ sphinx:
 formats: []
 
 python:
-  version: 3.7
+  version: 3.8
   install:
     - requirements: dev-requirements.txt
     - method: pip
diff --git a/dev-environment.yml b/dev-environment.yml
index bc0e45a1..cf19a330 100644
--- a/dev-environment.yml
+++ b/dev-environment.yml
@@ -2,7 +2,7 @@ name: xdem-dev
 channels:
   - conda-forge
 dependencies:
-  - python>=3.7
+  - python>=3.8
   - proj>=7.2
   - geopandas
   - numba
diff --git a/environment.yml b/environment.yml
index 121ddd5a..720c78d6 100644
--- a/environment.yml
+++ b/environment.yml
@@ -2,7 +2,7 @@ name: xdem
 channels:
   - conda-forge
 dependencies:
-  - python>=3.7
+  - python>=3.8
   - proj>=7.2
   - geopandas
   - numba
diff --git a/setup.py b/setup.py
index 922be0af..a0f17abe 100644
--- a/setup.py
+++ b/setup.py
@@ -39,7 +39,7 @@
         "opencv": ["opencv"],
         "pytransform3d": ["pytransform3d"],
     },
-    python_requires=">=3.7",
+    python_requires=">=3.8",
     scripts=[],
     zip_safe=False,
 )
diff --git a/xdem/terrain.py b/xdem/terrain.py
index 1dd4ca63..88d4bc8f 100644
--- a/xdem/terrain.py
+++ b/xdem/terrain.py
@@ -59,7 +59,7 @@ def _get_quadric_coefficients(
     L = resolution
 
     # Allocate the output.
-    output = np.empty((12,) + dem.shape, dtype=dem.dtype) + np.nan
+    output = np.full((12,) + dem.shape, fill_value=np.nan)
 
     # Convert the string to a number (fewer bytes to compare each iteration)
     if fill_method == "median":
@@ -337,7 +337,7 @@ def _get_windowed_indexes(
     """
 
     # Allocate the outputs.
-    output = np.empty((5,) + dem.shape, dtype=dem.dtype) + np.nan
+    output = np.full((5,) + dem.shape, fill_value=np.nan)
 
     # Half window size
     hw = int(np.floor(window_size / 2))
diff --git a/xdem/volume.py b/xdem/volume.py
index 895d24ef..08fb8802 100644
--- a/xdem/volume.py
+++ b/xdem/volume.py
@@ -71,7 +71,7 @@ def hypsometric_binning(ddem: np.ndarray, ref_dem: np.ndarray, bins: Union[float
     # Calculate statistics for each bin.
     # If no values exist, all stats should be nans (except count with should be 0)
     # medians, means, stds, nmads = (np.zeros(shape=bins.shape[0] - 1, dtype=ddem.dtype) * np.nan, ) * 4
-    values = np.zeros(shape=zbins.shape[0] - 1, dtype=ddem.dtype) * np.nan
+    values = np.full(shape=zbins.shape[0] - 1, fill_value=np.nan, dtype=ddem.dtype)
     counts = np.zeros_like(values, dtype=int)
     for i in np.arange(indices.min(), indices.max() + 1):
         values_in_bin = ddem[indices == i]
@@ -549,8 +549,8 @@ def get_regional_hypsometric_signal(ddem: Union[np.ndarray, np.ma.masked_array],
     unique_indices = np.unique(glacier_index_map)
 
     # Create empty (ddem) value and (pixel) count arrays which will be filled iteratively.
-    values = np.empty((n_bins, unique_indices.shape[0]), dtype=float) * np.nan
-    counts = np.empty((n_bins, unique_indices.shape[0]), dtype=float) * np.nan
+    values = np.full((n_bins, unique_indices.shape[0]), fill_value=np.nan, dtype=float)
+    counts = np.full((n_bins, unique_indices.shape[0]), fill_value=np.nan, dtype=float)
 
     # Start a counter of glaciers that are actually processed.
     count = 0