Commit 6a19fcba authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Refactored the term 'srf' to 'rsr'.


Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent e73589e3
......@@ -27,20 +27,20 @@ def RSR_reader(satellite, sensor, no_thermal=False, no_pan=False, v=False):
:param v: verbose mode
"""
# LayerBandsAssignment = META.get_LayerBandsAssignment(GMS_id, no_thermal=no_thermal, no_pan=no_pan)
SRF_dict = collections.OrderedDict()
SRF_dir = os.path.join(__path__[0], 'data', satellite, sensor)
RSR_dict = collections.OrderedDict()
RSR_dir = os.path.join(__path__[0], 'data', satellite, sensor)
for bandname in glob(os.path.join(SRF_dir, 'band_*')):
SRF_path = os.path.join(SRF_dir, bandname)
for bandname in glob(os.path.join(RSR_dir, 'band_*')):
RSR_path = os.path.join(RSR_dir, bandname)
try:
SRF_dict[bandname] = np.loadtxt(SRF_path, skiprows=1)
RSR_dict[bandname] = np.loadtxt(RSR_path, skiprows=1)
if v:
print('Reading SRF for %s %s, %s...' % (satellite, sensor, bandname))
print('Reading RSR for %s %s, %s...' % (satellite, sensor, bandname))
except FileNotFoundError:
raise FileNotFoundError('No spectral response functions found for %s %s %s at %s! >None< is returned.'
% (satellite, sensor, bandname, SRF_path))
% (satellite, sensor, bandname, RSR_path))
return SRF_dict
return RSR_dict
class RelativeSpectralResponse(object):
......@@ -65,9 +65,9 @@ class RelativeSpectralResponse(object):
if wvl_unit not in ['micrometers', 'nanometers']:
raise ValueError('Unknown wavelength unit %s.' % wvl_unit)
self.srfs_wvl = [] # wavelength positions with 1 nm precision
self.srfs = {} # srf values with 1 nm precision
self.srfs_norm01 = {} # srf values with 1 nm precision
self.rsr_wvl = [] # wavelength positions with 1 nm precision
self.rsrs = {} # RSR values with 1 nm precision
self.rsrs_norm01 = {} # RSR values with 1 nm precision
self.bands = []
self.wvl = []
self.wvl_unit = wvl_unit
......@@ -84,64 +84,64 @@ class RelativeSpectralResponse(object):
self.from_satellite_sensor(satellite, sensor)
def from_satellite_sensor(self, satellite, sensor):
srf_dict = RSR_reader(satellite, sensor, no_thermal=self.no_thermal, no_pan=self.no_pan,
rsr_dict = RSR_reader(satellite, sensor, no_thermal=self.no_thermal, no_pan=self.no_pan,
v=self.v) # type: collections.OrderedDict # (ordered according to LBA)
return self.from_dict(srf_dict)
return self.from_dict(rsr_dict)
def from_dict(self, srf_dict):
def from_dict(self, rsr_dict):
# type: (collections.OrderedDict) -> RelativeSpectralResponse
"""Create an instance of SRF from a dictionary.
"""Create an instance of RelativeSpectralResponse from a dictionary.
:param srf_dict: {'key_LayerBandsAssignment': <2D array: cols=[wvl,resp],rows=samples>}
:param rsr_dict: {'key_LayerBandsAssignment': <2D array: cols=[wvl,resp],rows=samples>}
"""
is_nm = [300 < np.max(srf_dict[band][:, 0]) < 15000 for band in srf_dict]
assert len(set(is_nm)) == 1, "'srf_dict' must contain only one wavelength unit."
is_nm = [300 < np.max(rsr_dict[band][:, 0]) < 15000 for band in rsr_dict]
assert len(set(is_nm)) == 1, "'rsr_dict' must contain only one wavelength unit."
scale_factor = 1 if is_nm[0] else 1000
all_wvls = np.concatenate(list((arr[:, 0] for key, arr in srf_dict.items()))) * scale_factor
all_wvls = np.concatenate(list((arr[:, 0] for key, arr in rsr_dict.items()))) * scale_factor
wvl = np.arange(all_wvls.min(), all_wvls.max() + self.specres_nm, self.specres_nm).astype(np.int16)
df = DataFrame(index=wvl)
bandnames = []
for band in srf_dict: # = OrderedDict -> order follows LayerBandsAssignment
for band in rsr_dict: # = OrderedDict -> order follows LayerBandsAssignment
bandname = band if not self.format_bandnames else ('B%s' % band if len(band) == 2 else 'B0%s' % band)
bandnames.append(bandname)
srfs = srf_dict[band][:, 1]
wvls = np.array(srf_dict[band][:, 0] * scale_factor)
rsrs = rsr_dict[band][:, 1]
wvls = np.array(rsr_dict[band][:, 0] * scale_factor)
# interpolate input SRFs to target spectral resolution
# interpolate input RSRs to target spectral resolution
wvl_s = wvl[np.abs(wvl - min(wvls)).argmin()]
wvl_e = wvl[np.abs(wvl - max(wvls)).argmin()]
wvls_tgtres = np.arange(wvl_s, wvl_e + self.specres_nm, self.specres_nm).astype(np.int16)
srfs_tgtres = interp1d(wvls, srfs, bounds_error=False, fill_value=0, kind='linear')(wvls_tgtres)
rsrs_tgtres = interp1d(wvls, rsrs, bounds_error=False, fill_value=0, kind='linear')(wvls_tgtres)
# join SRF for the current band to df
df = df.join(Series(srfs_tgtres, index=wvls_tgtres, name=bandname))
# join RSR for the current band to df
df = df.join(Series(rsrs_tgtres, index=wvls_tgtres, name=bandname))
df = df.fillna(0)
wvl = np.array(df.index.astype(np.float))
srfs_norm01 = df.to_dict(orient='series') # type: Dict[Series]
rsrs_norm01 = df.to_dict(orient='series') # type: Dict[Series]
################
# set attributes
################
for bN in srfs_norm01:
self.srfs_norm01[bN] = srf = np.array(srfs_norm01[bN], dtype=np.float)
self.srfs[bN] = srf / np.trapz(x=wvl, y=srf) # TODO seems like we NEED nanometers here; BUT WHY??
for bN in rsrs_norm01:
self.rsrs_norm01[bN] = rsr = np.array(rsrs_norm01[bN], dtype=np.float)
self.rsrs[bN] = rsr / np.trapz(x=wvl, y=rsr) # TODO seems like we NEED nanometers here; BUT WHY??
self.srfs_wvl = np.array(wvl)
self.rsr_wvl = np.array(wvl)
self.bands = bandnames
# FIXME this is not the GMS algorithm to calculate center wavelengths
# calculate center wavelengths
# TODO seems like we NEED nanometers here; BUT WHY??:
self.wvl = np.array([np.trapz(x=self.srfs_wvl, y=self.srfs_wvl * self.srfs[band]) for band in self.bands])
self.wvl = np.array([np.trapz(x=self.rsr_wvl, y=self.rsr_wvl * self.rsrs[band]) for band in self.bands])
# self.wvl = self.wvl if self.wvl_unit == 'micrometers' else np.array([int(i) for i in self.wvl])
self.srfs_wvl = self.srfs_wvl if self.wvl_unit == 'nanometers' else self.srfs_wvl / 1000
self.rsr_wvl = self.rsr_wvl if self.wvl_unit == 'nanometers' else self.rsr_wvl / 1000
self.wvl = self.wvl if self.wvl_unit == 'nanometers' else self.wvl / 1000
self.conv.update({key: value for key, value in zip(self.bands, self.wvl)})
......@@ -152,7 +152,7 @@ class RelativeSpectralResponse(object):
def instrument(self, bands):
return {
'rspf': np.vstack([self[band] for band in bands]),
'wvl_rsp': np.copy(self.srfs_wvl),
'wvl_rsp': np.copy(self.rsr_wvl),
'wvl_inst': np.copy(self.wvl),
'sol_irr': None
}
......@@ -160,26 +160,26 @@ class RelativeSpectralResponse(object):
def convert_wvl_unit(self):
"""Convert the wavelength unit to nanometers if they are in micrometers or vice versa."""
factor = 1/1000 if self.wvl_unit == 'nanometers' else 1000
self.srfs_wvl = self.srfs_wvl * factor
self.rsr_wvl = self.rsr_wvl * factor
self.wvl = self.wvl * factor
self.wvl_unit = 'nanometers' if self.wvl_unit == 'micrometers' else 'micrometers'
def __call__(self, band):
return self.srfs[band]
return self.rsrs[band]
def __getitem__(self, band):
return self.srfs[band]
return self.rsrs[band]
def __iter__(self):
for band in self.bands:
yield self[band]
def plot_srfs(self, figsize: tuple = (15, 5), band: Union[str, List[str]] = None, normalize: bool = True):
def plot_rsrs(self, figsize: tuple = (15, 5), band: Union[str, List[str]] = None, normalize: bool = True):
"""Show a plot of all spectral response functions.
:param figsize: figure size of the plot
:param band: band key to plot a single band instead of all bands
:param normalize: normalize SRFs to 0-1 (default: True)
:param normalize: normalize RSRs to 0-1 (default: True)
"""
if band and band not in self.bands:
raise ValueError("Parameter 'band' must be a string out of those: %s."
......@@ -188,8 +188,8 @@ class RelativeSpectralResponse(object):
plt.figure(figsize=figsize)
bands2plot = [band] if band else self.bands
for band in bands2plot:
srfs = list(self.srfs_norm01[band]) if normalize else list(self.srfs[band])
plt.plot(self.srfs_wvl, srfs, '-', label='Band %s' % band)
rsrs = list(self.rsrs_norm01[band]) if normalize else list(self.rsrs[band])
plt.plot(self.rsr_wvl, rsrs, '-', label='Band %s' % band)
plt.title(' '.join(['Spectral response functions', self.satellite, self.sensor]))
plt.xlabel('wavelength [%s]' % self.wvl_unit)
plt.ylabel('spectral response')
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment