Coverage for enpt/processors/atmospheric_correction/atmospheric_correction.py: 77%
137 statements
« prev ^ index » next coverage.py v7.4.1, created at 2024-03-07 11:39 +0000
1# -*- coding: utf-8 -*-
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/>.
30"""EnPT atmospheric correction module.
32Performs the atmospheric correction of EnMAP L1B data.
33"""
34import pprint
35import numpy as np
36from multiprocessing import cpu_count
37from logging import Logger
38import sys
40from packaging.version import parse as parse_version
41from sicor.sicor_enmap import sicor_ac_enmap
42from sicor.options import get_options as get_ac_options
44from ...model.images import EnMAPL1Product_SensorGeo
45from ...options.config import EnPTConfig
46from ...utils.path_generator import get_path_ac_options
48__author__ = 'Daniel Scheffler'
51class AtmosphericCorrector(object):
52 """Class for performing atmospheric correction of EnMAP L1 images using SICOR."""
54 def __init__(self, config: EnPTConfig = None):
55 """Create an instance of AtmosphericCorrector."""
56 self.cfg = config
58 def _get_sicor_options(self, enmap_ImageL1: EnMAPL1Product_SensorGeo, land_only=False) -> dict:
59 """Get a dictionary containing the SICOR options.
61 :param enmap_ImageL1: EnMAP Level 1 product in sensor geometry
62 :param land_only: True: SICOR is applied to land only; False: SICOR is applied to all pixels
63 :return: dictionary of SICOR options
64 """
65 path_opts = get_path_ac_options()
67 try:
68 options = get_ac_options(path_opts, validation=True)
70 # adjust options
71 if enmap_ImageL1.meta.aot is not None:
72 options["retrieval"]["default_aot_value"] = enmap_ImageL1.meta.aot
74 options["retrieval"]["cpu"] = self.cfg.CPUs or cpu_count()
75 options["retrieval"]["disable_progressbars"] = self.cfg.disable_progress_bars
77 # TODO: issue is closed -> revise
78 # temporarily disable uncertainty measures to avoid https://git.gfz-potsdam.de/EnMAP/sicor/-/issues/86
79 # if set to False, uncertainty values are not contained in the additional output of SICOR
80 options["retrieval"]["inversion"]["full"] = False
82 # set land_only mode
83 options["retrieval"]["land_only"] = land_only
85 # disable first guess water vapor retrieval for now
86 options["retrieval"]["state_vector"]["water_vapor"]["use_prior_mean"] = True
87 options["retrieval"]["state_vector"]["water_vapor"]["prior_mean"] = \
88 enmap_ImageL1.meta.water_vapour # = default = 2.5
90 # disable first guess liquid water retrieval for now
91 options["retrieval"]["state_vector"]["liquid_water"]["use_prior_mean"] = True
93 # disable first guess ice retrieval for now
94 options["retrieval"]["state_vector"]["ice"]["use_prior_mean"] = True
96 except FileNotFoundError:
97 raise FileNotFoundError(f'Could not locate options file for atmospheric correction at {path_opts}')
99 enmap_ImageL1.logger.debug('SICOR AC configuration: \n' +
100 pprint.pformat(options))
102 return options
104 def _is_acwater_operable(self, logger: Logger):
105 """Return True if ACWater/Polymer is operable, else raise a warning and return False."""
106 try:
107 import acwater as _acwater # noqa: F401
109 if parse_version(_acwater.__version__) < parse_version('0.3.0'):
110 if self.cfg.mode_ac in ['water', 'combined']:
111 logger.warning(f"The installed version of ACwater (v{_acwater.__version__}) is too old. "
112 f"At least version 0.3.0 is required. Instead of ACwater, SICOR is applied to water "
113 f"surfaces as a workaround.")
115 return False
117 except ImportError as e:
118 if self.cfg.mode_ac in ['water', 'combined']:
119 logger.warning(f"The atmospheric correction mode was set to '{self.cfg.mode_ac}' but "
120 f"ACwater cannot be imported (Error was: {e.msg}). "
121 f"As a fallback, SICOR is applied to water surfaces instead.")
122 return False
124 try:
125 if self.cfg.polymer_root:
126 sys.path.append(self.cfg.polymer_root)
127 from acwater.acwater import polymer_ac_enmap as _polymer_ac_enmap
128 if not _polymer_ac_enmap:
129 logger.warning("Polymer is not callable. "
130 "As a fallback, SICOR is applied to water surfaces instead.")
131 return False
132 except ImportError as e:
133 if self.cfg.mode_ac in ['water', 'combined']:
134 logger.warning(f"The atmospheric correction mode was set to '{self.cfg.mode_ac}' but "
135 f"Polymer cannot be imported (Error was: {e.msg}). "
136 f"As a fallback, SICOR is applied to water surfaces instead.")
137 return False
139 return True
141 def _run_ac__land_mode(self,
142 enmap_ImageL1: EnMAPL1Product_SensorGeo
143 ) -> (np.ndarray, np.ndarray, dict):
144 """Run atmospheric correction in 'land' mode, i.e., use SICOR for all surfaces."""
145 if 2 in enmap_ImageL1.vnir.mask_landwater[:]:
146 enmap_ImageL1.logger.info(
147 "Running atmospheric correction in 'land' mode, i.e., SICOR is applied to ALL surfaces. "
148 "Uncertainty is expected for water surfaces because SICOR is designed for land only.")
150 # run SICOR
151 # NOTE: - boa_ref_vnir, boa_ref_swir: reflectance between 0 and 1
152 # - res: a dictionary containing retrieval maps with path lengths of the three water phases
153 # and several retrieval uncertainty measures
154 # -> cwv_model, liq_model, ice_model, intercept_model, slope_model, toa_model,
155 # sx (posterior predictive uncertainty matrix), scem (correlation error matrix),
156 # srem (relative error matrix)
157 # optional (if options["retrieval"]["inversion"]["full"] = True):
158 # jacobian, convergence, iterations, gain, averaging_kernel, cost_function,
159 # dof (degrees of freedom), information_content, retrieval_noise, smoothing_error
160 # -> SWIR geometry (?) # FIXME
161 boa_ref_vnir, boa_ref_swir, land_additional_results = \
162 sicor_ac_enmap(enmap_l1b=enmap_ImageL1,
163 options=self._get_sicor_options(enmap_ImageL1, land_only=False),
164 unknowns=True,
165 logger=enmap_ImageL1.logger)
167 return boa_ref_vnir, boa_ref_swir, land_additional_results
169 def _run_ac__water_mode(self, enmap_ImageL1: EnMAPL1Product_SensorGeo
170 ) -> (np.ndarray, np.ndarray):
171 """Run atmospheric correction in 'water' mode, i.e., use ACWater/Polymer for water surfaces only.
173 NOTE:
174 - Land surfaces are NOT included in the EnMAP L2A product.
175 - The output radiometric unit for water surfaces is 'water leaving reflectance'.
176 """
177 from acwater.acwater import polymer_ac_enmap
179 if 1 in enmap_ImageL1.vnir.mask_landwater[:]:
180 enmap_ImageL1.logger.info(
181 "Running atmospheric correction in 'water' mode, i.e., ACWater/Polymer is applied to water "
182 "surfaces only. Note that land surfaces will NOT be included in the EnMAP L2A product.")
184 # run ACWater/Polymer for water surfaces only
185 # NOTE: polymer_ac_enmap() returns masked (nan) values for land
186 # - res: a dictionary containing retrieval maps with several additional retrieval measures
187 # -> chla, bitmask, logfb, Rnir, Rgli
188 try:
189 wl_ref_vnir, wl_ref_swir, water_additional_results = \
190 polymer_ac_enmap(enmap_l1b=enmap_ImageL1,
191 config=self.cfg,
192 detector='vnir')
194 # Overwrite SWIR with 0 for water pixels (POLYMER does not produce a SWIR output)
195 # and NaNs for all other pixels (NaNs are later set to no-data)
196 # -> not needed anymore if implemented in ACwater - https://gitlab.awi.de/phytooptics/acwater/-/issues/23
197 wl_ref_swir = np.zeros_like(wl_ref_swir)
198 wl_ref_swir[enmap_ImageL1.swir.mask_landwater[:] != 2] = np.nan
200 except: # noqa
201 enmap_ImageL1.logger.error(
202 "The atmospheric correction for water surfaces based on ACwater/Polymer failed (issue tracker at "
203 "https://gitlab.awi.de/phytooptics/acwater/-/issues).\n"
204 "Alternatively, you may run EnPT in the 'land' atmospheric correction mode based on SICOR.\n"
205 "The error message is now raised:"
206 )
207 raise
209 return wl_ref_vnir, wl_ref_swir, water_additional_results
211 def _run_ac__combined_mode(self,
212 enmap_ImageL1: EnMAPL1Product_SensorGeo
213 ) -> (np.ndarray, np.ndarray, dict):
214 """Run atmospheric corr. in 'combined' mode, i.e., use SICOR for land and ACWater/Polymer for water surfaces.
216 NOTE:
217 - The output radiometric units are:
218 - 'surface reflectance' for land surfaces
219 - 'water leaving reflectance' for water surfaces
220 - There might be visible edge effects, e.g., at coastlines.
221 """
222 from acwater.acwater import polymer_ac_enmap
224 only = 'water' if 1 not in enmap_ImageL1.vnir.mask_landwater[:] else 'land'
225 if 1 not in enmap_ImageL1.vnir.mask_landwater[:] or \
226 2 not in enmap_ImageL1.vnir.mask_landwater[:]:
227 enmap_ImageL1.logger.info(
228 f"Running atmospheric correction in 'combined' mode, i.e., SICOR is applied to land and "
229 f"ACWater/Polymer is applied to water surfaces. But the input image only contains {only} surfaces.")
231 # run SICOR for land surfaces only
232 boa_ref_vnir_land, boa_ref_swir_land, land_additional_results = \
233 sicor_ac_enmap(enmap_l1b=enmap_ImageL1,
234 options=self._get_sicor_options(enmap_ImageL1, land_only=True),
235 unknowns=True,
236 logger=enmap_ImageL1.logger)
238 # run ACWater/Polymer for water surfaces only
239 # NOTE: polymer_ac_enmap() returns masked (nan) values for land
240 # - res: a dictionary containing retrieval maps with several additional retrieval measures
241 # -> chla, bitmask, logfb, Rnir, Rgli
242 try:
243 wl_ref_vnir_water, wl_ref_swir_water, water_additional_results = \
244 polymer_ac_enmap(enmap_l1b=enmap_ImageL1,
245 config=self.cfg,
246 detector='vnir')
248 # Overwrite SWIR with 0 for water pixels (POLYMER does not produce a SWIR output)
249 # and NaNs for all other pixels (NaNs are later set to no-data)
250 # -> not needed anymore if implemented in ACwater - https://gitlab.awi.de/phytooptics/acwater/-/issues/23
251 wl_ref_swir_water = np.zeros_like(wl_ref_swir_water)
252 wl_ref_swir_water[enmap_ImageL1.swir.mask_landwater[:] != 2] = np.nan
254 except: # noqa
255 enmap_ImageL1.logger.error(
256 "The atmospheric correction for water surfaces based on ACwater/Polymer failed (issue tracker at "
257 "https://gitlab.awi.de/phytooptics/acwater/-/issues).\n"
258 "Alternatively, you may run EnPT in the 'land' atmospheric correction mode based on SICOR.\n"
259 "The error message is now raised:"
260 )
261 raise
263 # Overwrite the SICOR output at water positions with the output from ACwater/Polymer
264 # -> output contains Water-leaving-reflectance over water and BOA-reflectance over land
265 water_mask_vnir = enmap_ImageL1.vnir.mask_landwater[:] == 2 # 2 = water
266 water_mask_swir = enmap_ImageL1.swir.mask_landwater[:] == 2 # 2 = water
267 wlboa_ref_vnir = np.where(water_mask_vnir[:, :, None], wl_ref_vnir_water, boa_ref_vnir_land)
268 wlboa_ref_swir = np.where(water_mask_swir[:, :, None], wl_ref_swir_water, boa_ref_swir_land)
270 return wlboa_ref_vnir, wlboa_ref_swir, water_additional_results, land_additional_results
272 @staticmethod
273 def _validate_ac_results(reflectance_vnir: np.ndarray,
274 reflectance_swir: np.ndarray,
275 logger: Logger):
276 for detectordata, detectorname in zip([reflectance_vnir, reflectance_swir],
277 ['VNIR', 'SWIR']):
278 mean0 = np.nanmean(detectordata[:, :, 0])
279 std0 = np.nanstd(detectordata[:, :, 0])
281 if np.isnan(mean0) or \
282 mean0 == 0 or \
283 std0 == 0:
284 logger.warning(f'The atmospheric correction returned empty {detectorname} bands!')
286 def run_ac(self,
287 enmap_ImageL1: EnMAPL1Product_SensorGeo
288 ) -> EnMAPL1Product_SensorGeo:
289 """Run atmospheric correction according to the specified 'mode_ac' parameter.
291 :param enmap_ImageL1: input EnMAP image containing TOA reflectance (an instance EnMAPL1Product_SensorGeo)
292 :return: atmospherically corrected output EnMAP image containing BOA reflectance / water leaving reflectance
293 (an instance EnMAPL1Product_SensorGeo)
294 """
295 enmap_ImageL1.set_SWIRattr_with_transformedVNIRattr('mask_landwater', src_nodata=0, tgt_nodata=0)
297 enmap_ImageL1.logger.info(
298 f"Starting atmospheric correction for VNIR and SWIR detector in '{self.cfg.mode_ac}' mode. "
299 f"Source radiometric unit code is '{enmap_ImageL1.meta.vnir.unitcode}'.")
301 # set initial value water_additional_results
302 water_additional_results = None
304 # run the AC
305 acwater_operable = self._is_acwater_operable(enmap_ImageL1.logger)
307 if self.cfg.mode_ac in ['water', 'combined'] and not acwater_operable:
308 # use SICOR as fallback AC for water surfaces if ACWater/Polymer is not installed
309 reflectance_vnir, reflectance_swir, land_additional_results = \
310 self._run_ac__land_mode(enmap_ImageL1)
312 else:
313 if self.cfg.mode_ac == 'combined':
314 reflectance_vnir, reflectance_swir, water_additional_results, land_additional_results = \
315 self._run_ac__combined_mode(enmap_ImageL1)
317 elif self.cfg.mode_ac == 'water':
318 reflectance_vnir, reflectance_swir, water_additional_results = \
319 self._run_ac__water_mode(enmap_ImageL1)
321 elif self.cfg.mode_ac == 'land':
322 reflectance_vnir, reflectance_swir, land_additional_results = \
323 self._run_ac__land_mode(enmap_ImageL1)
325 else:
326 raise ValueError(self.cfg.mode_ac,
327 "Unexpected 'mode_ac' parameter. "
328 "Choose one out of 'land', 'water', 'combined'.")
330 # validate outputs
331 self._validate_ac_results(reflectance_vnir, reflectance_swir, logger=enmap_ImageL1.logger)
333 # join results
334 enmap_ImageL1.logger.info('Joining results of atmospheric correction.')
336 for in_detector, out_detector in zip([enmap_ImageL1.vnir, enmap_ImageL1.swir],
337 [reflectance_vnir, reflectance_swir]):
338 data_ac_scaled_float = out_detector * self.cfg.scale_factor_boa_ref
340 if self.cfg.mode_ac in ['combined', 'water'] and acwater_operable:
341 # Overwrite the POLYMER output at land positions (set to NaN by POLYMER)
342 # AND other pixels which may also be set to NaN by POLYMER with the output nodata value
343 data_ac_scaled_float = np.nan_to_num(data_ac_scaled_float, nan=self.cfg.output_nodata_value)
345 # NOTE: - geotransform and projection are missing due to sensor geometry
346 # - remaining NaNs not due to POLYMER intentionally cause a numpy warning when casting to int16
347 in_detector.data = data_ac_scaled_float.astype(np.int16)
349 in_detector.detector_meta.unit = '0-%d' % self.cfg.scale_factor_boa_ref
350 in_detector.detector_meta.unitcode = 'BOARef'
352 # FIXME: Consider to also join SICOR's land_additional_results
353 # (contains three phases of water maps and several retrieval uncertainty measures)
355 # join additional results from ACwater/Polymer
356 if water_additional_results and self.cfg.polymer_additional_results:
358 water_mask = enmap_ImageL1.vnir.mask_landwater[:] == 2
359 for k in water_additional_results.keys():
360 if k == 'polymer_bitmask':
361 # the bitmask already explicitly indicates land pixels with "1"
362 continue
363 else:
364 v = water_additional_results[k]
365 v[~water_mask] = -9999
366 v[np.isnan(v)] = -9999
368 water_additional_results[k] = v
370 enmap_ImageL1.vnir.polymer_logchl = water_additional_results['polymer_logchl']
371 enmap_ImageL1.vnir.polymer_logfb = water_additional_results['polymer_logfb']
372 enmap_ImageL1.vnir.polymer_rgli = water_additional_results['polymer_rgli']
373 enmap_ImageL1.vnir.polymer_rnir = water_additional_results['polymer_rnir']
374 enmap_ImageL1.vnir.polymer_bitmask = water_additional_results['polymer_bitmask']
376 return enmap_ImageL1