From c5cf192c31c84171055b353f79d1fd3dce2317f8 Mon Sep 17 00:00:00 2001
From: Daniel Scheffler <danschef@gfz-potsdam.de>
Date: Tue, 5 Mar 2024 13:06:49 +0100
Subject: [PATCH] Replace interp1d in metadata_sensorgeo.py, reader.py, and
 spatial_transform.py.

Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>
---
 enpt/io/reader.py                                 |  9 ++++-----
 enpt/model/metadata/metadata_sensorgeo.py         | 13 ++++++-------
 .../spatial_transform/spatial_transform.py        | 15 ++++++++-------
 3 files changed, 18 insertions(+), 19 deletions(-)

diff --git a/enpt/io/reader.py b/enpt/io/reader.py
index de0432d5..732bd8bb 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 a95da215..d51ddaf4 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/processors/spatial_transform/spatial_transform.py b/enpt/processors/spatial_transform/spatial_transform.py
index 17dcbf05..3ec62842 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
 
-- 
GitLab