Coverage for enpt/model/srf.py: 71%

1# -*- coding: utf-8 -*- 

2 

3# EnPT, EnMAP Processing Tool - A Python package for pre-processing of EnMAP Level-1B data 

4# 

5# Copyright (C) 2018-2024 Karl Segl (GFZ Potsdam, segl@gfz-potsdam.de), Daniel Scheffler 

6# (GFZ Potsdam, danschef@gfz-potsdam.de), Niklas Bohn (GFZ Potsdam, nbohn@gfz-potsdam.de), 

7# Stéphane Guillaso (GFZ Potsdam, stephane.guillaso@gfz-potsdam.de) 

8# 

9# This software was developed within the context of the EnMAP project supported 

10# by the DLR Space Administration with funds of the German Federal Ministry of 

11# Economic Affairs and Energy (on the basis of a decision by the German Bundestag: 

12# 50 EE 1529) and contributions from DLR, GFZ and OHB System AG. 

13# 

14# This program is free software: you can redistribute it and/or modify it under 

15# the terms of the GNU General Public License as published by the Free Software 

16# Foundation, either version 3 of the License, or (at your option) any later 

17# version. Please note the following exception: `EnPT` depends on tqdm, which 

18# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files 

19# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore". 

20# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE. 

21# 

22# This program is distributed in the hope that it will be useful, but WITHOUT 

23# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 

24# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more 

25# details. 

26# 

27# You should have received a copy of the GNU Lesser General Public License along 

28# with this program. If not, see <https://www.gnu.org/licenses/>. 

29 

30"""EnPT module for handling spectral response functions.""" 

31 

32from typing import Union, List 

33 

34import numpy as np 

35from matplotlib import pyplot as plt 

36from scipy import stats 

37 

38__author__ = 'Daniel Scheffler' 

39 

40 

41class SRF(object): 

42 def __init__(self, wvl_unit: str = 'nanometers', wvl_min: float = 400, wvl_max: float = 2500, specres_nm: float = 1, 

43 format_bandnames: bool = False, v: bool = False): 

44 """SRF instance provides SRF functions, wavelength positions, etc.. 

45 

46 :param wvl_unit: the wavelengths unit to be used within SRF instance ('nanometers' or 'micrometers) 

47 :param wvl_min: 

48 :param wvl_max: 

49 :param specres_nm: output spectral resolution of SRFs in nanometers 

50 :param format_bandnames: whether to format default strings from LayerBandsAssignment as 'B01', 'B02' etc.. 

51 :param v: verbose mode 

52 """ 

53 if wvl_unit not in ['micrometers', 'nanometers']: 

54 raise ValueError('Unknown wavelength unit %s.' % wvl_unit) 

55 

56 self.srfs_wvl = [] # wavelength positions with 1 nm precision 

57 self.srfs = {} # srf values with 1 nm precision 

58 self.srfs_norm01 = {} # srf values with 1 nm precision 

59 self.bands = [] 

60 self.wvl = None 

61 self.wvl_unit = wvl_unit 

62 self.wvl_min = wvl_min 

63 self.wvl_max = wvl_max 

64 self.specres_nm = specres_nm 

65 self.format_bandnames = format_bandnames 

66 self.conv = {} 

67 self.v = v 

68 

69 @staticmethod 

70 def compute_gaussian_srf(cwl: float, fwhm: float, wvl_min: float, wvl_max: float, wvl_res: float, 

71 normalize: bool = True) -> np.ndarray: 

72 """Compute a spectral response function based on center wavelength and bandwidth using on a gaussian curve. 

73 

74 :param cwl: target center wavelength position 

75 :param fwhm: target bandwidth (full width half maximum) 

76 :param wvl_min: minimum wavelength to compute spectral response for 

77 :param wvl_max: maximum wavelength to compute spectral response for 

78 :param wvl_res: spectral resolution at which spectral response is to be computed 

79 :param normalize: whether to normalize the output spectral response to values between 0 and 1 

80 :return: 2D numpy.ndarray: rows: response per wavelength; columns: wavelength/response 

81 """ 

82 x = np.arange(wvl_min, wvl_max, wvl_res) 

83 dist = stats.norm(cwl, fwhm) 

84 

85 with np.errstate(under='ignore'): 

86 # This suppresses Numpy warnings: "underflow encountered in exp/divide/multiply" due to very small values 

87 y = dist.pdf(x) 

88 

89 if normalize: 

90 y *= (1.0 / y.max()) 

91 

92 rsp = np.empty((x.size, 2), dtype=float) 

93 rsp[:, 0] = x 

94 rsp[:, 1] = y 

95 

96 return rsp 

97 

98 @classmethod 

99 def from_cwl_fwhm(cls, cwls: Union[list, np.ndarray], fwhms: Union[list, np.ndarray], **kwargs: dict) -> 'SRF': 

100 """Create an instance of SRF based on center wavelength positions and bandwidths (using gaussian responses). 

101 

102 :param cwls: center wavelength positions 

103 :param fwhms: bandwidths 

104 :param kwargs: Keyword arguments to be passed to SRF.__init__(). 

105 :return: SRF instance 

106 """ 

107 srf = cls(**kwargs) 

108 

109 if srf.wvl_unit != 'nanometers': 

110 cwls, fwhms = cwls * 1000, fwhms * 1000 

111 

112 srf.bands = [str(i + 1) for i in range(len(cwls))] 

113 

114 # compute SRFs and set attributes 

115 for bN, cwl, fwhm in zip(srf.bands, cwls, fwhms): 

116 gaussian_srf = cls.compute_gaussian_srf(cwl, fwhm, srf.wvl_min, srf.wvl_max, srf.specres_nm) 

117 srf.srfs_wvl = gaussian_srf[:, 0].flatten() 

118 srf_norm01 = gaussian_srf[:, 1].flatten() 

119 srf.srfs_norm01[bN] = srf_norm01 

120 with np.errstate(under='ignore'): # suppress Warning: "underflow encountered in divide due to small values 

121 srf.srfs[bN] = srf_norm01 / np.trapz(x=srf.srfs_wvl, y=srf_norm01) 

122 

123 srf.wvl = np.array(cwls) 

124 

125 srf.conv.update({key: value for key, value in zip(srf.bands, srf.wvl)}) 

126 srf.conv.update({value: key for key, value in zip(srf.bands, srf.wvl)}) 

127 

128 return srf 

129 

130 def instrument(self, bands): 

131 return { 

132 'rspf': np.vstack([self[band] for band in bands]), 

133 'wvl_rsp': np.copy(self.srfs_wvl), 

134 'wvl_inst': np.copy(self.wvl), 

135 'sol_irr': None 

136 } 

137 

138 def convert_wvl_unit(self): 

139 """Convert the wavelength unit to nanometers if they are in micrometers or vice versa.""" 

140 factor = 1/1000 if self.wvl_unit == 'nanometers' else 1000 

141 self.srfs_wvl = self.srfs_wvl * factor 

142 self.wvl = self.wvl * factor 

143 self.wvl_unit = 'nanometers' if self.wvl_unit == 'micrometers' else 'micrometers' 

144 

145 def __call__(self, band): 

146 return self.srfs[band] 

147 

148 def __getitem__(self, band): 

149 return self.srfs[band] 

150 

151 def __iter__(self): 

152 for band in self.bands: 

153 yield self[band] 

154 

155 def plot_srfs(self, figsize: tuple = (15, 5), band: Union[str, List[str]] = None, normalize: bool = True): 

156 """Show a plot of all spectral response functions. 

157 

158 :param figsize: figure size of the plot 

159 :param band: band key to plot a single band instead of all bands 

160 :param normalize: normalize SRFs to 0-1 (default: True) 

161 """ 

162 if band and band not in self.bands: 

163 raise ValueError("Parameter 'band' must be a string out of those: %s." 

164 % ', '.join(self.bands)) 

165 

166 plt.figure(figsize=figsize) 

167 bands2plot = [band] if band else self.bands 

168 for band in bands2plot: 

169 srfs = list(self.srfs_norm01[band]) if normalize else list(self.srfs[band]) 

170 plt.plot(self.srfs_wvl, srfs, '-', label='Band %s' % band) 

171 plt.title('EnMAP spectral response functions') 

172 plt.xlabel('wavelength [%s]' % self.wvl_unit) 

173 plt.ylabel('spectral response') 

174 plt.legend(loc='upper right')