Coverage for spechomo/resampling.py : 87%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
3# spechomo, Spectral homogenization of multispectral satellite data
4#
5# Copyright (C) 2019-2021
6# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
7# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,
8# Germany (https://www.gfz-potsdam.de/)
9#
10# This software was developed within the context of the GeoMultiSens project funded
11# by the German Federal Ministry of Education and Research
12# (project grant code: 01 IS 14 010 A-C).
13#
14# Licensed under the Apache License, Version 2.0 (the "License");
15# you may not use this file except in compliance with the License.
16# You may obtain a copy of the License at
17#
18# http://www.apache.org/licenses/LICENSE-2.0
19#
20# Please note the following exception: `spechomo` depends on tqdm, which is
21# distributed under the Mozilla Public Licence (MPL) v2.0 except for the files
22# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore".
23# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE.
24#
25# Unless required by applicable law or agreed to in writing, software
26# distributed under the License is distributed on an "AS IS" BASIS,
27# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
28# See the License for the specific language governing permissions and
29# limitations under the License.
31from multiprocessing import Pool
32from typing import Union, TYPE_CHECKING # noqa F401 # flake8 issue
34import numpy as np
35from geoarray import GeoArray
36from geoarray.baseclasses import get_array_tilebounds
38from .logging import SpecHomo_Logger
40if TYPE_CHECKING:
41 from pyrsr import RSR # noqa F401 # flake8 issue
44class SpectralResampler(object):
45 """Class for spectral resampling of a single spectral signature (1D-array) or an image (3D-array)."""
47 def __init__(self, wvl_src, rsr_tgt, logger=None):
48 # type: (np.ndarray, RSR, str) -> None
49 """Get an instance of the SpectralResampler class.
51 :param wvl_src: center wavelength positions of the source spectrum
52 :param rsr_tgt: relative spectral response of the target instrument as an instance of
53 pyrsr.RelativeSpectralRespnse.
54 """
55 # privates
56 self._wvl_1nm = None
57 self._rsr_1nm = {}
59 wvl_src = np.array(wvl_src, dtype=float).flatten()
60 if rsr_tgt.wvl_unit != 'nanometers':
61 rsr_tgt.convert_wvl_unit()
63 self.wvl_src_nm = wvl_src if max(wvl_src) > 100 else wvl_src * 1000
64 self.rsr_tgt = rsr_tgt
65 self.logger = logger or SpecHomo_Logger(__name__) # must be pickable
67 @property
68 def wvl_1nm(self):
69 # spectral resampling of input image to 1 nm resolution
70 if self._wvl_1nm is None:
71 self._wvl_1nm = np.arange(np.ceil(self.wvl_src_nm.min()), np.floor(self.wvl_src_nm.max()), 1)
72 return self._wvl_1nm
74 @property
75 def rsr_1nm(self):
76 if not self._rsr_1nm:
77 from scipy.interpolate import interp1d # import here to avoid static TLS ImportError
79 for band in self.rsr_tgt.bands:
80 # resample RSR to 1 nm
81 self._rsr_1nm[band] = \
82 interp1d(self.rsr_tgt.rsrs_wvl, self.rsr_tgt.rsrs[band],
83 bounds_error=False, fill_value=0, kind='linear')(self.wvl_1nm)
85 # validate
86 assert len(self._rsr_1nm[band]) == len(self.wvl_1nm)
88 return self._rsr_1nm
90 def resample_signature(self, spectrum, scale_factor=10000, nodataVal=None, alg_nodata='radical', v=False):
91 # type: (np.ndarray, int, Union[int, float], str, bool) -> np.ndarray
92 """Resample the given spectrum according to the spectral response functions of the target instument.
94 :param spectrum: spectral signature data
95 :param scale_factor: the scale factor to apply to the given spectrum when it is plotted (default: 10000)
96 :param nodataVal: no data value to be respected during resampling
97 :param alg_nodata: algorithm how to deal with pixels where the spectral bands of the source image
98 contain nodata within the spectral response of a target band
99 'radical': set output band to nodata
100 'conservative': use existing spectral information and ignore nodata
101 (might alter the outpur spectral information,
102 e.g., at spectral absorption bands)
103 :param v: enable verbose mode (shows a plot of the resampled spectrum) (default: False)
104 :return: resampled spectral signature
105 """
106 if not spectrum.ndim == 1:
107 raise ValueError("The array of the given spectral signature must be 1-dimensional. "
108 "Received a %s-dimensional array." % spectrum.ndim)
109 spectrum = np.array(spectrum, dtype=float).flatten()
110 assert spectrum.size == self.wvl_src_nm.size
112 # convert spectrum to one multispectral image pixel and resample it
113 spectrum_rsp = \
114 self.resample_image(spectrum.reshape(1, 1, spectrum.size),
115 tiledims=(1, 1),
116 nodataVal=nodataVal,
117 alg_nodata=alg_nodata,
118 CPUs=1)\
119 .ravel()
121 if v:
122 from matplotlib import pyplot as plt
124 plt.figure()
125 for band in self.rsr_tgt.bands:
126 plt.plot(self.wvl_1nm, self.rsr_1nm[band] / max(self.rsr_1nm[band]))
127 plt.plot(self.wvl_src_nm, spectrum / scale_factor, '.')
128 plt.plot(self.rsr_tgt.wvl, spectrum_rsp / scale_factor, 'x', color='r')
129 plt.show()
131 return spectrum_rsp
133 def resample_spectra(self, spectra, chunksize=200, nodataVal=None, alg_nodata='radical', CPUs=None):
134 # type: (Union[GeoArray, np.ndarray], int, Union[int, float], str, Union[None, int]) -> np.ndarray
135 """Resample the given spectral signatures according to the spectral response functions of the target instrument.
137 :param spectra: spectral signatures, provided as 2D array
138 (rows: spectral samples, columns: spectral information / bands)
139 :param chunksize: defines how many spectral signatures are resampled per CPU
140 :param nodataVal: no data value to be respected during resampling
141 :param alg_nodata: algorithm how to deal with pixels where the spectral bands of the source image
142 contain nodata within the spectral response of a target band
143 'radical': set output band to nodata
144 'conservative': use existing spectral information and ignore nodata
145 (might alter the outpur spectral information,
146 e.g., at spectral absorption bands)
147 :param CPUs: CPUs to use for processing
148 """
149 # input validation
150 if not spectra.ndim == 2:
151 ValueError("The given spectra array must be 2-dimensional. Received a %s-dimensional array."
152 % spectra.ndim)
153 assert spectra.shape[1] == self.wvl_src_nm.size
155 # convert spectra to one multispectral image column
156 im_col = spectra.reshape(spectra.shape[0], 1, spectra.shape[1])
158 im_col_rsp = self.resample_image(im_col,
159 tiledims=(1, chunksize),
160 nodataVal=nodataVal,
161 alg_nodata=alg_nodata,
162 CPUs=CPUs)
163 spectra_rsp = im_col_rsp.reshape(im_col_rsp.shape[0], im_col_rsp.shape[2])
165 return spectra_rsp
167 def resample_image(self, image_cube, tiledims=(20, 20), nodataVal=None, alg_nodata='radical', CPUs=None):
168 # type: (Union[GeoArray, np.ndarray], tuple, Union[int, float], str, Union[None, int]) -> np.ndarray
169 """Resample the given spectral image cube according to the spectral response functions of the target instrument.
171 :param image_cube: image (3D array) containing the spectral information in the third dimension
172 :param tiledims: dimension of tiles to be used during computation (rows, columns)
173 :param nodataVal: nodata value of the input image
174 :param alg_nodata: algorithm how to deal with pixels where the spectral bands of the source image
175 contain nodata within the spectral response of a target band
176 'radical': set output band to nodata
177 'conservative': use existing spectral information and ignore nodata
178 (might alter the output spectral information,
179 e.g., at spectral absorption bands)
180 :param CPUs: CPUs to use for processing
181 :return: resampled spectral image cube
182 """
183 # input validation
184 if not image_cube.ndim == 3:
185 raise ValueError("The given image cube must be 3-dimensional. Received a %s-dimensional array."
186 % image_cube.ndim)
187 assert image_cube.shape[2] == self.wvl_src_nm.size
189 image_cube = GeoArray(image_cube)
191 (R, C), B = image_cube.shape[:2], len(self.rsr_tgt.bands)
192 image_rsp = np.zeros((R, C, B), dtype=image_cube.dtype)
194 initargs = image_cube, self.rsr_tgt, self.rsr_1nm, self.wvl_src_nm, self.wvl_1nm
195 if CPUs is None or CPUs > 1:
196 with Pool(CPUs, initializer=_initializer_mp, initargs=initargs) as pool:
197 tiles_rsp = pool.starmap(_resample_tile_mp,
198 [(bounds, nodataVal, alg_nodata)
199 for bounds in get_array_tilebounds(image_cube.shape, tiledims)])
200 pool.close() # needed to make coverage work in multiprocessing
201 pool.join()
203 else:
204 _initializer_mp(*initargs)
205 tiles_rsp = [_resample_tile_mp(bounds, nodataVal, alg_nodata)
206 for bounds in get_array_tilebounds(image_cube.shape, tiledims)]
208 for ((rS, rE), (cS, cE)), tile_rsp in tiles_rsp:
209 image_rsp[rS: rE + 1, cS: cE + 1, :] = tile_rsp
211 return image_rsp
214_globs_mp = dict()
217def _initializer_mp(image_cube, rsr_tgt, rsr_1nm, wvl_src_nm, wvl_1nm):
218 # type: (np.ndarray, RSR, RSR, np.ndarray, np.ndarray) -> None
219 global _globs_mp
220 _globs_mp.update(dict(
221 image_cube=image_cube,
222 rsr_tgt=rsr_tgt,
223 rsr_1nm=rsr_1nm,
224 wvl_src_nm=wvl_src_nm,
225 wvl_1nm=wvl_1nm,
226 ))
229def _resample_tile_mp(tilebounds, nodataVal=None, alg_nodata='radical'):
230 # TODO speed up by using numba as described here https://krstn.eu/fast-linear-1D-interpolation-with-numba/
232 from scipy.interpolate import interp1d
234 (rS, rE), (cS, cE) = tilebounds
236 # get global share variables
237 tiledata = _globs_mp['image_cube'][rS: rE + 1, cS: cE + 1, :] # type: Union[GeoArray, np.ndarray]
238 rsr_tgt = _globs_mp['rsr_tgt']
239 rsr_1nm = _globs_mp['rsr_1nm']
240 wvl_src_nm = _globs_mp['wvl_src_nm']
241 wvl_1nm = _globs_mp['wvl_1nm']
243 if alg_nodata not in ['radical', 'conservative']:
244 raise ValueError(alg_nodata)
246 tile_rsp = np.zeros((*tiledata.shape[:2], len(rsr_tgt.bands)), dtype=tiledata.dtype)
248 if nodataVal is not None:
249 if np.isfinite(nodataVal):
250 _mask_bool3d = tiledata == nodataVal
251 mask_anynodata = np.any(_mask_bool3d, axis=2)
252 mask_allnodata = np.all(_mask_bool3d, axis=2)
253 else:
254 raise NotImplementedError(nodataVal)
256 tiledata = tiledata.astype(float)
257 tiledata[_mask_bool3d] = np.nan
259 # fill pixels with all bands nodata
260 tile_rsp[mask_allnodata] = nodataVal
262 # upsample spectra containing data to 1nm and get nan-mask
263 tilespectra_1nm = interp1d(wvl_src_nm, tiledata[~mask_allnodata],
264 axis=1, bounds_error=False, fill_value=0, kind='linear')(wvl_1nm)
265 isnan = np.isnan(tilespectra_1nm)
266 nan_in_spec = np.any(isnan, axis=1)
268 # compute resampled values for pixels without nodata
269 tilespectra_1nm_nonan = tilespectra_1nm[~nan_in_spec, :]
270 for band_idx, band in enumerate(rsr_tgt.bands):
271 # compute the resampled image cube (np.average computes the weighted mean value)
272 if rsr_1nm[band].sum() != 0:
273 res = np.average(tilespectra_1nm_nonan, weights=rsr_1nm[band], axis=1)
274 # NOTE: rounding here is important to prevent -9999. converted to -9998
275 if not np.issubdtype(tile_rsp.dtype, np.floating):
276 res = np.around(res)
277 else:
278 res = nodataVal
280 tile_rsp[~mask_anynodata, band_idx] = res
282 # compute resampled values for pixels with nodata
283 tilespectra_1nm_withnan = tilespectra_1nm[nan_in_spec, :]
284 tilespectra_1nm_withnan_ma = np.ma.masked_invalid(tilespectra_1nm_withnan)
285 mask_withnan = mask_anynodata & ~mask_allnodata
286 isnan_sub = isnan[nan_in_spec, :]
288 for band_idx, band in enumerate(rsr_tgt.bands):
289 res_ma = np.ma.average(tilespectra_1nm_withnan_ma, axis=1, weights=rsr_1nm[band])
291 # in case all 1nm bands for the current band are NaN, the resulting average will also be NaN => fill
292 res = np.ma.filled(res_ma, nodataVal)
294 if alg_nodata == 'radical':
295 # set those output values to nodata where the input bands within the target RSR contain any nodata
296 badspec = np.any(isnan_sub & (rsr_1nm[band] > 0), axis=1)
297 res[badspec] = nodataVal
299 if not np.issubdtype(tile_rsp.dtype, np.floating):
300 res = np.around(res)
302 tile_rsp[mask_withnan, band_idx] = res
304 else:
305 # spectral resampling of input image to 1 nm resolution
306 tile_1nm = interp1d(wvl_src_nm, tiledata,
307 axis=2, bounds_error=False, fill_value=0, kind='linear')(wvl_1nm)
309 for band_idx, band in enumerate(rsr_tgt.bands):
310 # compute the resampled image cube (np.average computes the weighted mean value)
311 res = np.average(tile_1nm, weights=rsr_1nm[band], axis=2)
312 # NOTE: rounding here is important to prevent -9999. converted to -9998
313 tile_rsp[:, :, band_idx] = res if np.issubdtype(tile_rsp.dtype, np.floating) else np.around(res)
315 return tilebounds, tile_rsp