diff --git a/HISTORY.rst b/HISTORY.rst
index cb4f1998c554f969dfc3f7c11e917f5a81c23ed4..8ed4aacce0807d108a3a410e98704277c274864b 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,6 +2,12 @@
History
=======
+0.20.1 (2024-03-05)
+-------------------
+
+* !105: Replaced legacy scipy.interpolate.interp1d with similar 1-dimensional interpolation functions.
+
+
0.20.0 (2024-02-09)
-------------------
diff --git a/enpt/cli.py b/enpt/cli.py
index cc4f30cf17361d9c0c2c39938f53cf3c381f94f9..bd8336cbddd659aeece5f78678712c2f63582885 100644
--- a/enpt/cli.py
+++ b/enpt/cli.py
@@ -135,7 +135,7 @@ def get_enpt_argparser():
help="Algorithm for dead pixel correction ('spectral' or 'spatial')")
add('--deadpix_P_interp_spectral', type=str, default="linear",
help="Spectral interpolation algorithm to be used during dead pixel correction "
- "('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic')")
+ "('linear', 'quadratic', 'cubic')")
add('--deadpix_P_interp_spatial', type=str, default="linear",
help="Spatial interpolation algorithm to be used during dead pixel correction "
"('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic')")
diff --git a/enpt/io/reader.py b/enpt/io/reader.py
index de0432d5d4703ab3f9aa5619f907e780e880e445..732bd8bbb6c67e090640543192aa4615af506f62 100644
--- a/enpt/io/reader.py
+++ b/enpt/io/reader.py
@@ -33,7 +33,6 @@ import logging
import os
from fnmatch import filter
import numpy as np
-from scipy.interpolate import interp1d
from ..model.images import EnMAPL1Product_SensorGeo
from ..options.config import EnPTConfig
@@ -162,14 +161,14 @@ class L1B_Reader(object):
pass
-def Solar_Irradiance_reader(path_solar_irr_model: str, resol_nm: float = None, wvl_min_nm: float = None,
- wvl_max_nm: float = None) -> np.ndarray:
+def read_solar_irradiance(path_solar_irr_model: str, resol_nm: float = None, wvl_min_nm: float = None,
+ wvl_max_nm: float = None) -> np.ndarray:
"""Read the given solar irradiance file and return an array of irradiances.
:param path_solar_irr_model: file path to solar irradiance model
-> must be arranged like that:
- col0 = Wavelength[nm]; col1 = Solar Irradiance [W/m2/µm])
+ col0 = Wavelength[nm]; col1 = Solar Irradiance [W/m2/µm]
:param resol_nm: spectral resolution for returned irradiances [nanometers]
:param wvl_min_nm: minimum wavelength of returned irradiances [nanometers]
:param wvl_max_nm: maximum wavelength of returned irradiances [nanometers]
@@ -180,5 +179,5 @@ def Solar_Irradiance_reader(path_solar_irr_model: str, resol_nm: float = None, w
wvl_min = (np.min(sol_irr[:, 0]) if wvl_min_nm is None else wvl_min_nm)
wvl_max = (np.max(sol_irr[:, 0]) if wvl_max_nm is None else wvl_max_nm)
wvl_rsp = np.arange(wvl_min, wvl_max, resol_nm)
- sol_irr = interp1d(sol_irr[:, 0], sol_irr[:, 1], kind='linear')(wvl_rsp)
+ sol_irr = np.interp(wvl_rsp, sol_irr[:, 0], sol_irr[:, 1])
return sol_irr
diff --git a/enpt/model/metadata/metadata_sensorgeo.py b/enpt/model/metadata/metadata_sensorgeo.py
index a95da215ed2ac49fc320095a490ba4317cd56c79..d51ddaf450b182d8cdf07998823eab5535a5a39b 100644
--- a/enpt/model/metadata/metadata_sensorgeo.py
+++ b/enpt/model/metadata/metadata_sensorgeo.py
@@ -384,11 +384,10 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object):
return gA
else:
- from scipy.interpolate import interp1d
- coeffs_interp = \
- interp1d(gA.meta.band_meta['wavelength'], gA[:],
- axis=2, kind='linear', fill_value="extrapolate", bounds_error=False)(self.wvl_center)
-
+ from scipy.interpolate import make_interp_spline
+ # equivalent to legacy scipy.interpolate.interp1d
+ # interp1d(wvl, gA[:], axis=2, kind='linear', fill_value="extrapolate", bounds_error=False)(self.wvl_center)
+ coeffs_interp = make_interp_spline(gA.meta.band_meta['wavelength'], gA[:], k=1, axis=2)(self.wvl_center)
gA_ = GeoArray(coeffs_interp)
gA_.meta.band_meta['wavelength'] = self.wvl_center
@@ -444,11 +443,11 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object):
return lons, lats
def calc_solar_irradiance_CWL_FWHM_per_band(self) -> np.array:
- from ...io.reader import Solar_Irradiance_reader
+ from ...io.reader import read_solar_irradiance
self.logger.debug('Calculating solar irradiance...')
- sol_irr = Solar_Irradiance_reader(path_solar_irr_model=self.cfg.path_solar_irr, wvl_min_nm=350, wvl_max_nm=2500)
+ sol_irr = read_solar_irradiance(path_solar_irr_model=self.cfg.path_solar_irr, wvl_min_nm=350, wvl_max_nm=2500)
irr_bands = []
for band in self.srf.bands:
diff --git a/enpt/options/config.py b/enpt/options/config.py
index c5d391510f92688cee7d572ae6853baa448520bd..cbb8e0e1bc8e48bf5a0ea652affa77e092bb222e 100644
--- a/enpt/options/config.py
+++ b/enpt/options/config.py
@@ -322,7 +322,7 @@ class EnPTConfig(object):
:key deadpix_P_interp_spectral:
Spectral interpolation algorithm to be used during dead pixel correction
- ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic')
+ ('linear', 'quadratic', 'cubic')
:key deadpix_P_interp_spatial:
Spatial interpolation algorithm to be used during dead pixel correction
diff --git a/enpt/options/options_default.json b/enpt/options/options_default.json
index 1096e7c821841662cb6ed7b22c8702b79307de6b..c16a5022866ffbf22276f1c74fbc99b98659fc61 100644
--- a/enpt/options/options_default.json
+++ b/enpt/options/options_default.json
@@ -84,7 +84,7 @@
'spectral': interpolate in the spectral domain
'spatial': interpolate in the spatial domain*/
"interp_method_spectral": "linear", /*spectral interpolation algorithm
- ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic', ...)*/
+ ('linear', 'quadratic', 'cubic')*/
"interp_method_spatial": "linear" /*spatial interpolation algorithm
('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic', ...)*/
},
diff --git a/enpt/options/options_schema.py b/enpt/options/options_schema.py
index 3b832b35696d37a03a9f9dca7cd9a060a242a782..20bbb2749a6579327de0ff98da5d23deee454a53 100644
--- a/enpt/options/options_schema.py
+++ b/enpt/options/options_schema.py
@@ -106,7 +106,7 @@ enpt_schema_input = dict(
run_processor=dict(type='boolean', required=False),
algorithm=dict(type='string', required=False, allowed=['spectral', 'spatial']),
interp_method_spectral=dict(type='string', required=False,
- allowed=['linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic']),
+ allowed=['linear', 'quadratic', 'cubic']),
interp_method_spatial=dict(type='string', required=False,
allowed=['linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic']),
)),
diff --git a/enpt/processors/dead_pixel_correction/dead_pixel_correction.py b/enpt/processors/dead_pixel_correction/dead_pixel_correction.py
index 89f04df15624b0adc4119e4dac5863df0694ade7..69b28cc688c419a65d017b38c0adbdff29398a59 100644
--- a/enpt/processors/dead_pixel_correction/dead_pixel_correction.py
+++ b/enpt/processors/dead_pixel_correction/dead_pixel_correction.py
@@ -34,13 +34,15 @@ Performs the Dead Pixel Correction using a linear interpolation in spectral dime
from typing import Union
from numbers import Number # noqa: F401
import logging
+from warnings import filterwarnings
import numpy as np
import numpy_indexed as npi
from multiprocessing import Pool, cpu_count
-from scipy.interpolate import griddata, interp1d
-from pandas import DataFrame
-from geoarray import GeoArray
+from scipy.interpolate import griddata, make_interp_spline
+filterwarnings("ignore", "\nPyArrow", DeprecationWarning) # mute pandas warning
+from pandas import DataFrame # noqa: E402
+from geoarray import GeoArray # noqa: E402
__author__ = 'Daniel Scheffler'
@@ -70,8 +72,7 @@ class Dead_Pixel_Corrector(object):
:param algorithm: algorithm how to correct dead pixels
'spectral': interpolate in the spectral domain
'spatial': interpolate in the spatial domain
- :param interp_spectral: spectral interpolation algorithm
- (‘linear’, ‘nearest’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’, etc.)
+ :param interp_spectral: spectral interpolation algorithm (‘linear’, ‘quadratic’, ‘cubic’)
:param interp_spatial: spatial interpolation algorithm ('linear', 'bilinear', 'cubic', 'spline')
:param CPUs: number of CPUs to use for interpolation (only relevant if algorithm = 'spatial')
:param logger:
@@ -106,7 +107,7 @@ class Dead_Pixel_Corrector(object):
raise ValueError("Dead pixel map and image to be corrected must have equal shape.")
image_corrected = interp_nodata_along_axis(image2correct, axis=2, nodata=deadpixel_map[:],
- method=self.interp_alg_spectral, fill_value='extrapolate')
+ method=self.interp_alg_spectral)
return image_corrected
@@ -117,7 +118,7 @@ class Dead_Pixel_Corrector(object):
if deadpixel_map.shape != image2correct.shape:
raise ValueError("Dead pixel map and image to be corrected must have equal shape.")
- band_indices_with_nodata = np.argwhere(np.any(np.any(deadpixel_map, axis=0), axis=0)).flatten()
+ band_indices_with_nodata = np.argwhere(np.any(np.any(deadpixel_map[:], axis=0), axis=0)).flatten()
image_sub = image2correct[:, :, band_indices_with_nodata]
deadpixel_map_sub = deadpixel_map[:, :, band_indices_with_nodata]
@@ -136,7 +137,7 @@ class Dead_Pixel_Corrector(object):
# correct remaining nodata by spectral interpolation (e.g., outermost columns)
if np.isnan(image2correct).any():
image2correct = interp_nodata_along_axis(image2correct, axis=2, nodata=np.isnan(image2correct),
- method=self.interp_alg_spectral, fill_value='extrapolate')
+ method=self.interp_alg_spectral)
return image2correct
@@ -205,21 +206,24 @@ def interp_nodata_along_axis_2d(data_2d: np.ndarray,
axis: int = 0,
nodata: Union[np.ndarray, Number] = np.nan,
method: str = 'linear',
- fill_value: Union[float, str] = 'extrapolate'):
- """Interpolate a 2D array along the given axis (based on scipy.interpolate.interp1d).
+ **kw):
+ """Interpolate a 2D array along the given axis (based on scipy.interpolate.make_interp_spline).
:param data_2d: data to interpolate
:param axis: axis to interpolate (0: along columns; 1: along rows)
:param nodata: nodata array in the shape of data or nodata value
- :param method: interpolation method (‘linear’, ‘nearest’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’, etc.)
- :param fill_value: value to fill into positions where no interpolation is possible
- - if 'extrapolate': extrapolate the missing values
+ :param method: interpolation method (‘linear’, ‘quadratic’, ‘cubic’)
+ :param kw: keyword arguments to be passed to scipy.interpolate.make_interp_spline
:return: interpolated array
"""
if data_2d.ndim != 2:
raise ValueError('Expected a 2D array. Received a %dD array.' % data_2d.ndim)
if axis > data_2d.ndim:
raise ValueError("axis=%d is out of bounds for data with %d dimensions." % (axis, data_2d.ndim))
+ if method not in ['linear', 'quadratic', 'cubic']:
+ raise ValueError(f"'{method}' is not a valid interpolation method. "
+ f"Choose between 'linear', 'quadratic', and 'cubic'.")
+ degree = 1 if method == 'linear' else 2 if method == 'quadratic' else 3
data_2d = data_2d if axis == 1 else data_2d.T
@@ -240,10 +244,8 @@ def interp_nodata_along_axis_2d(data_2d: np.ndarray,
if goodpos.size > 1:
data_2d_grouped_rows = data_2d[indices_unique_rows]
-
data_2d_grouped_rows[:, badpos] = \
- interp1d(goodpos, data_2d_grouped_rows[:, goodpos],
- axis=1, kind=method, fill_value=fill_value, bounds_error=False)(badpos)
+ make_interp_spline(goodpos, data_2d_grouped_rows[:, goodpos], axis=1, k=degree, **kw)(badpos)
data_2d[indices_unique_rows, :] = data_2d_grouped_rows
@@ -254,15 +256,14 @@ def interp_nodata_along_axis(data,
axis=0,
nodata: Union[np.ndarray, Number] = np.nan,
method: str = 'linear',
- fill_value: Union[float, str] = 'extrapolate'):
- """Interpolate a 2D or 3D array along the given axis (based on scipy.interpolate.interp1d).
+ **kw):
+ """Interpolate a 2D or 3D array along the given axis (based on scipy.interpolate.make_interp_spline).
:param data: data to interpolate
:param axis: axis to interpolate (0: along columns; 1: along rows, 2: along bands)
:param nodata: nodata array in the shape of data or nodata value
- :param method: interpolation method (‘linear’, ‘nearest’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’, etc.)
- :param fill_value: value to fill into positions where no interpolation is possible
- - if 'extrapolate': extrapolate the missing values
+ :param method: interpolation method (‘linear’, 'quadratic', 'cubic')
+ :param kw: keyword arguments to be passed to scipy.interpolate.make_interp_spline
:return: interpolated array
"""
assert axis <= 2
@@ -272,7 +273,7 @@ def interp_nodata_along_axis(data,
raise ValueError('No-data mask and data must have the same shape.')
if data.ndim == 2:
- return interp_nodata_along_axis_2d(data, axis=axis, nodata=nodata, method=method, fill_value=fill_value)
+ return interp_nodata_along_axis_2d(data, axis=axis, nodata=nodata, method=method, **kw)
else:
def reshape_input(In):
@@ -293,7 +294,7 @@ def interp_nodata_along_axis(data,
data_2d=reshape_input(data),
nodata=reshape_input(nodata) if isinstance(nodata, np.ndarray) else nodata,
axis=axis if axis != 2 else 1,
- method=method, fill_value=fill_value))
+ method=method, **kw))
def interp_nodata_spatially_2d(data_2d: np.ndarray,
diff --git a/enpt/processors/spatial_transform/spatial_transform.py b/enpt/processors/spatial_transform/spatial_transform.py
index 17dcbf0581b2e446d37d92032c912075c663bd99..3ec62842f664df633d80aa4b527527e4139be72b 100644
--- a/enpt/processors/spatial_transform/spatial_transform.py
+++ b/enpt/processors/spatial_transform/spatial_transform.py
@@ -32,7 +32,7 @@ from typing import Union, Tuple, List, Optional, Sequence # noqa: F401
from multiprocessing import Pool, cpu_count
from collections import OrderedDict
import numpy as np
-from scipy.interpolate import interp1d, LinearNDInterpolator
+from scipy.interpolate import LinearNDInterpolator, make_interp_spline
from scipy.spatial import Delaunay
from geoarray import GeoArray
from natsort import natsorted
@@ -405,26 +405,27 @@ class RPC_Geolayer_Generator(object):
raise ValueError(arr.ndim, '2D numpy array expected.')
if along_axis not in [0, 1]:
raise ValueError(along_axis, "The 'axis' parameter must be set to 0 or 1")
-
- kw = dict(kind='linear', fill_value='extrapolate')
+ arr[0, 0] = np.nan
+ arr[0, 1] = np.nan
+ arr[1, 1] = np.nan
nans = np.isnan(arr)
if along_axis == 0:
- # extrapolate in top/bottom direction
+ # extrapolate linearly in top/bottom direction
cols_with_nan = np.arange(arr.shape[1])[np.any(nans, axis=0)]
for col in cols_with_nan:
data = arr[:, col]
idx_goodvals = np.argwhere(~nans[:, col]).flatten()
- arr[:, col] = interp1d(idx_goodvals, data[idx_goodvals], **kw)(range(data.size))
+ arr[:, col] = make_interp_spline(idx_goodvals, data[idx_goodvals], k=1)(range(data.size))
else:
- # extrapolate in left/right direction
+ # extrapolate linearly in left/right direction
rows_with_nan = np.arange(arr.shape[0])[np.any(nans, axis=1)]
for row in rows_with_nan:
data = arr[row, :]
idx_goodvals = np.argwhere(~nans[row, :]).flatten()
- arr[row, :] = interp1d(idx_goodvals, data[idx_goodvals], **kw)(range(data.size))
+ arr[row, :] = make_interp_spline(idx_goodvals, data[idx_goodvals], k=1)(range(data.size))
return arr
diff --git a/tests/test_dead_pixel_correction.py b/tests/test_dead_pixel_correction.py
index 75eb2691fd8b4f2625c1183013020e9224c1a8e0..cc233b14b18687ad546a78affcafebd035f960e2 100644
--- a/tests/test_dead_pixel_correction.py
+++ b/tests/test_dead_pixel_correction.py
@@ -126,10 +126,6 @@ class Test_interp_nodata_along_axis_2d(TestCase):
data_int = interp_nodata_along_axis_2d(self.get_data2d(), axis=0, nodata=mask_nodata, method='linear')
assert np.array_equal(data_int, arr_exp), 'Computed %s.' % data_int
- data_int = interp_nodata_along_axis_2d(self.get_data2d(), axis=0, method='linear', fill_value=-1)
- arr_exp = np.array([[0, 0, 2], [3, 5, 5], [-1, 10, 8]])
- assert np.array_equal(data_int, arr_exp), 'Computed %s.' % data_int
-
def test_axis_1(self):
data_int = interp_nodata_along_axis_2d(self.get_data2d(), axis=1, method='linear')
arr_exp = np.array([[0, 0, 2], [3, 4, 5], [12, 10, 8]])
@@ -139,15 +135,13 @@ class Test_interp_nodata_along_axis_2d(TestCase):
data_int = interp_nodata_along_axis_2d(self.get_data2d(), axis=1, nodata=mask_nodata, method='linear')
assert np.array_equal(data_int, arr_exp), 'Computed %s.' % data_int
- data_int = interp_nodata_along_axis_2d(self.get_data2d(), axis=1, method='linear', fill_value=-1)
- arr_exp = np.array([[0, 0, 2], [3, 4, 5], [-1, 10, 8]])
- assert np.array_equal(data_int, arr_exp), 'Computed %s.' % data_int
-
def test_bad_args(self):
with pytest.raises(ValueError):
interp_nodata_along_axis_2d(self.get_data2d(), axis=3)
with pytest.raises(ValueError):
interp_nodata_along_axis_2d(np.dstack([self.get_data2d(), self.get_data2d()]))
+ with pytest.raises(ValueError):
+ interp_nodata_along_axis_2d(self.get_data2d(), method='unsupported_method')
class Test_interp_nodata_along_axis(TestCase):