diff --git a/AUTHORS.rst b/AUTHORS.rst index 1fa5394515a51dda0b378da6a6c0d9f5196d0574..512816a75b29b15926790a665ad49fe626192e53 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -13,6 +13,7 @@ Contributors * Daniel Scheffler (main developer of the EnPT source code) * Niklas Bohn - (main developer of the SICOR atmopsheric correction source code) + (main developer of the SICOR atmospheric correction source code) * André Hollstein * Stéphane Guillaso +* Brenner Silva diff --git a/HISTORY.rst b/HISTORY.rst index 57d96ad014f7a8575e18458afb48d2dc435f57c7..d60a4f555954da8ed236cd6d6e6e5c5ebc69e9a2 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -19,6 +19,10 @@ History * Don't assert existing file paths in config validation n case IS_ENPT_GUI_TEST==1. * Default orthorecification algorithm is now 'gauss'. +* Added config parameters to run EnPT in 3 AC modes: 'land', 'water', 'combined'. +* Added some boilerplate code in atmospheric_correction.py which is to be replaced by separate AC calls for water and + land surfaces later. + 0.17.2 (2021-03-04) ------------------- diff --git a/README.rst b/README.rst index 0e975b4d606b0ba582a5e323883e2bf7c46f1fa1..351582e34feab708861eb1263d08427513ee6cd4 100644 --- a/README.rst +++ b/README.rst @@ -31,6 +31,7 @@ Feature overview * radiometric conversion to top-of-atmosphere radiance * dead pixel correction * atmospheric correction (based on SICOR_) +* atmospheric correction of water surfaces (based on Polymer_ via `ACwater Polymer`_). * conversion of top-of-atmosphere-radiance to top-of-atmosphere-reflectance * detection and correction of geometric misregistrations compared to user provided spatial reference (based on AROSICS_) * orthorectification @@ -93,4 +94,6 @@ This package was created with Cookiecutter_ and the `audreyr/cookiecutter-pypack .. _coverage: https://enmap.git-pages.gfz-potsdam.de/GFZ_Tools_EnMAP_BOX/EnPT/coverage/ .. _nosetests: https://enmap.git-pages.gfz-potsdam.de/GFZ_Tools_EnMAP_BOX/EnPT/nosetests_reports/nosetests.html .. _SICOR: https://git.gfz-potsdam.de/EnMAP/sicor -.. _AROSICS: https://git.gfz-potsdam.de/danschef/arosics \ No newline at end of file +.. _AROSICS: https://git.gfz-potsdam.de/danschef/arosics +.. _`ACwater Polymer`: https://gitlab.awi.de/phytooptics/acwater +.. _Polymer: https://forum.hygeos.com diff --git a/docs/about.rst b/docs/about.rst index 80aa9c935fc4f4f7838b0356d9814d708d085a5e..055a2b75eea548c67286c27af72efeb0a08ec271 100644 --- a/docs/about.rst +++ b/docs/about.rst @@ -16,10 +16,12 @@ Feature overview * read EnMAP Level-1B input data * radiometric conversion to top-of-atmosphere radiance * dead pixel correction -* atmospheric correction (based on SICOR_) +* atmospheric correction (based on SICOR_ for land and `ACwater Polymer`_ via Polymer_ for water surfaces) * detection and correction of geometric misregistrations compared to user provided spatial reference (based on AROSICS_) * orthorectification * write EnMAP Level-2 output data .. _SICOR: https://git.gfz-potsdam.de/EnMAP/sicor .. _AROSICS: https://git.gfz-potsdam.de/danschef/arosics +.. _`ACwater Polymer`: https://gitlab.awi.de/phytooptics/acwater +.. _Polymer: https://forum.hygeos.com diff --git a/docs/algorithm_descriptions.rst b/docs/algorithm_descriptions.rst index 9f9a1fdbfb475e75b500aa3924435306e437f5fb..e9368e1be372c93bda855ae29abb6d19c7db111e 100644 --- a/docs/algorithm_descriptions.rst +++ b/docs/algorithm_descriptions.rst @@ -99,8 +99,17 @@ to BOA- (bottom-of-atmosphere / surface) reflectance. SICOR is a Python based op German Research Centre for Geosciences (GFZ) Potsdam. For details on the underlying algorithm, please refer to the `documentation pages of SICOR`_. +Optionally, EnPT retrieves water reflectance above the surface using `ACwater Polymer`_. +ACwater Polymer is a "wrapper" package (developed at the Alfred-Wegener-Institute, Bremerhaven) +for the `Polymer`_ atmospheric correction (AC) algorithm (developed by Hygeos, Inc). +Polymer AC is based on an optimization technique that considers atmospheric and oceanic signals to retrieve +normalized spectral reflectance above water. For details regarding the Polymer algorithm, +users are referred to `Steinmetz F, Deschamps P-Y, Ramon R., Opt. Express. 2011; 19`__. +__ https://doi.org/10.1364/OE.19.009783 +.. _`ACwater Polymer`: https://gitlab.awi.de/phytooptics/acwater +.. _Polymer: https://www.hygeos.com/polymer Spatial Co-Registration *********************** diff --git a/docs/installation.rst b/docs/installation.rst index 29196c5579cac2858b8131e68e8532d056a20316..baecc6a5a4abcf54ec11343003b1dc33a5e121ca 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -70,6 +70,35 @@ you through the process. +Optional: Install ACwater for advanced atmospheric correction over water surfaces +--------------------------------------------------------------------------------- + +For atmospheric correction of water surfaces using the Polymer algorithm instead of SICOR_ (which is mainly +designed for land surfaces), the additional package ACwater_ (a Polymer wrapper developed by AWI Bremerhaven) +is required. + +1. Using a previously created enpt conda environment (as described above), first install some dependencies: + + .. code-block:: bash + + $ conda activate enpt + $ conda install -c conda-forge cdsapi cython gdal netcdf4 pygrib pyhdf xarray + $ pip install ecmwf-api-client + +2. Then register at the `HYGEOS support forum`_, download polymer_ from there, unpack it and + run :code:`pip install .` from the root directory of Polymer. + +3. Finally install ACwater: + + .. code-block:: bash + + $ pip install git+https://gitlab.awi.de/phytooptics/acwater.git + + +Further details about the installation of ACwater can be found in the `ACwater Polymer installation instructions`_. + + + .. note:: EnPT has been tested with Python 3.6+., i.e., should be fully compatible to all Python versions from 3.6 onwards. @@ -78,3 +107,8 @@ you through the process. .. _pip: https://pip.pypa.io .. _Python installation guide: http://docs.python-guide.org/en/latest/starting/installation/ .. _conda: https://conda.io/docs +.. _ACwater: https://gitlab.awi.de/phytooptics/acwater/ +.. _`ACwater Polymer installation instructions`: https://gitlab.awi.de/phytooptics/acwater/-/blob/master/docs/installation.rst +.. _HYGEOS support forum: https://forum.hygeos.com +.. _`polymer`: https://forum.hygeos.com +.. _SICOR: https://git.gfz-potsdam.de/EnMAP/sicor diff --git a/enpt/cli.py b/enpt/cli.py index 5a6cffc57a6ec7981d7d321d8aed9f26b925e3f0..a2e673cd1aeefe424c7d319db77727e377848277 100644 --- a/enpt/cli.py +++ b/enpt/cli.py @@ -106,13 +106,27 @@ def get_enpt_argparser(): help='Enable VNIR/SWIR co-registration') add('--path_reference_image', type=str, default=None, help='Reference image for co-registration.') + add('--polymer_root', type=str, default=None, + help='Polymer root directory (that contains the subdirectory for ancillary data)') add('--enable_ac', type=_str2bool, default=True, nargs='?', const=True, help="Enable atmospheric correction using SICOR algorithm (default: True). If False, the L2A output contains " "top-of-atmosphere reflectance") + add('--mode_ac', type=str, default=None, nargs='?', + help="3 modes to determine which atmospheric correction is applied at which surfaces (default: land): " + "('land', water', 'combined')") + add('--auto_download_ecmwf', type=_str2bool, default=True, nargs='?', const=True, + help='Automatically download ECMWF data for atmospheric correction') + add('--enable_ice_retrieval', type=_str2bool, default=True, nargs='?', const=True, + help='Enable ice retrieval (default); increases accuracy of water vapour retrieval') add('--enable_cloud_screening', type=_str2bool, default=False, nargs='?', const=True, help='Enable cloud screening during atmospheric correction') add('--scale_factor_boa_ref', type=int, default=10000, help='Scale factor to be applied to BOA reflectance result') + add('--threads', type=int, default=-1, + help='Number of threads in ACwater Polymer: 0 for single thread; < 0 for as many as there are CPUs; ' + 'and > 0 gives the number of threads') + add('--blocksize', type=int, default=100, + help='Block size in ACwater Polymer') add('--run_smile_P', type=_str2bool, default=False, nargs='?', const=True, help='Enable extra smile detection and correction (provider smile coefficients are ignored)') add('--run_deadpix_P', type=_str2bool, default=True, nargs='?', const=True, diff --git a/enpt/model/images/images_mapgeo.py b/enpt/model/images/images_mapgeo.py index d015dd9e6de8da8f7ed5c9cad81b8be5c2d8177e..9610a60d3355074dbcd92b3cf165b662d17d37d9 100644 --- a/enpt/model/images/images_mapgeo.py +++ b/enpt/model/images/images_mapgeo.py @@ -143,7 +143,7 @@ class EnMAPL2Product_MapGeo(_EnMAP_Image): @property def logger(self) -> EnPT_Logger: - """Get a an instance of enpt.utils.logging.EnPT_Logger. + """Get an instance of enpt.utils.logging.EnPT_Logger. NOTE: - The logging level will be set according to the user inputs of EnPT. diff --git a/enpt/options/config.py b/enpt/options/config.py index 6ab6a6145526f33bae242a7d30b60d81175cae8f..ac9249dc393cac0da69fb314744196dd185dd72a 100644 --- a/enpt/options/config.py +++ b/enpt/options/config.py @@ -62,6 +62,56 @@ __author__ = 'Daniel Scheffler' path_enptlib = os.path.dirname(pkgutil.get_loader("enpt").path) path_options_default = os.path.join(path_enptlib, 'options', 'options_default.json') +try: + # from acwater.acwater import polymer_ac_enmap + path_polymer = os.path.abspath(os.path.join(os.path.dirname(pkgutil.get_loader("polymer").path), os.pardir)) +except AttributeError: + path_polymer = '' + +config_for_testing_water = dict( + path_l1b_enmap_image=os.path.abspath( + os.path.join(path_enptlib, '..', 'tests', 'data', 'EnMAP_Level_1B', + # Arcachon + 'ENMAP01-____L1B-DT000400126_20170218T110115Z_002_V000204_20200206T182719Z__rows700-730.zip' + + # Arcachon full tile 2 + # 'ENMAP01-____L1B-DT000400126_20170218T110115Z_002_V000204_20200206T182719Z.zip' + )), + # path_l1b_enmap_image_gapfill=os.path.abspath( + # os.path.join(path_enptlib, '..', 'tests', 'data', 'EnMAP_Level_1B', + # 'ENMAP01-____L1B-DT000400126_20170218T110115Z_002_V000204_20200206T182719Z__rows700-730.zip')), + path_dem=os.path.abspath( + os.path.join(path_enptlib, '..', 'tests', 'data', + 'ENMAP01-____L1B-DT000400126_20170218T110115Z_002_V000204_20200206T182719Z__tile2' + '__DEM_ASTER.bsq')), + log_level='DEBUG', + output_dir=os.path.join(path_enptlib, '..', 'tests', 'data', 'test_outputs'), + disable_progress_bars=False, + is_dummy_dataformat=False, + auto_download_ecmwf=True, + average_elevation=0, + deadpix_P_algorithm='spectral', + deadpix_P_interp_spatial='linear', + deadpix_P_interp_spectral='linear', + enable_cloud_screening=False, + enable_ice_retrieval=True, + enable_keystone_correction=False, + enable_vnir_swir_coreg=False, + n_lines_to_append=None, + ortho_resampAlg='gauss', + run_deadpix_P=True, + run_smile_P=False, + scale_factor_boa_ref=10000, + scale_factor_toa_ref=10000, + enable_ac=True, + mode_ac='combined', + polymer_root=path_polymer, + threads=-1, + blocksize=100, + vswir_overlap_algorithm='swir_only', + CPUs=16 +) + config_for_testing = dict( path_l1b_enmap_image=os.path.abspath( @@ -135,6 +185,8 @@ config_for_testing_dlr = dict( enable_absolute_coreg=True, path_reference_image=os.path.join(path_enptlib, '..', 'tests', 'data', 'T30TXQ_20170218T110111_B05__sub.tif'), enable_ac=True, + mode_ac='land', + enable_ice_retrieval=False, CPUs=32, ortho_resampAlg='gauss', vswir_overlap_algorithm='swir_only' @@ -218,16 +270,42 @@ class EnPTConfig(object): :key path_reference_image: Reference image for co-registration. + :key polymer root: + Polymer root directory (that contains the subdirectory for ancillary data). + :key enable_ac: Enable atmospheric correction using SICOR algorithm (default: True). If False, the L2A output contains top-of-atmosphere reflectance. + :key mode_ac: + 3 modes to determine which atmospheric correction is applied at which surfaces (default: land): + - 'land': SICOR (developed for land surfaces is applied to land AND water surfaces + - 'water': POLYMER (developed for water surfaces) is applied to water only + (land surfaces are no included in the L2A product) + - 'combined': SICOR is applied to land and POLYMER is applied to water surfaces; + NOTE that this may result in edge effects, e.g., at coastlines + + :key auto_download_ecmwf: + Automatically download ECMWF data for atmospheric correction + + :key enable_ice_retrieval: + Enable ice retrieval (default); increases accuracy of water vapour retrieval + :key enable_cloud_screening: Enable cloud screening during atmospheric correction :key scale_factor_boa_ref: Scale factor to be applied to BOA reflectance result + :key threads: + number of threads for multiprocessing of blocks (see bellow): + - 'threads = 0': for single thread + - 'threads < 0': for as many threads as there are CPUs + - 'threads > 0': gives the number of threads + + :key blocksize: + block size for multiprocessing + :key run_smile_P: Enable extra smile detection and correction (provider smile coefficients are ignored) @@ -321,9 +399,15 @@ class EnPTConfig(object): self.path_reference_image = gp('path_reference_image') # atmospheric_correction + self.polymer_root = gp('polymer_root') self.enable_ac = gp('enable_ac') + self.mode_ac = gp('mode_ac') + self.auto_download_ecmwf = gp('auto_download_ecmwf') + self.enable_ice_retrieval = gp('enable_ice_retrieval') self.enable_cloud_screening = gp('enable_cloud_screening') self.scale_factor_boa_ref = gp('scale_factor_boa_ref') + self.threads = gp('threads') + self.blocksize = gp('blocksize') # smile self.run_smile_P = gp('run_smile_P') diff --git a/enpt/options/options_default.json b/enpt/options/options_default.json index 2356ce568e6e7a8701d08fc21cde042fd48ce705..8311e1ec73492f5e48716aafb57c8c66ae5aae22 100644 --- a/enpt/options/options_default.json +++ b/enpt/options/options_default.json @@ -52,10 +52,26 @@ }, "atmospheric_correction": { + "polymer_root": "", /*Polymer root directory (that contains the subdirectory for ancillary data)*/ "enable_ac": true, /*enable atmospheric correction using SICOR algorithm (default: True). If False, the L2A output contains top-of-atmosphere reflectance.*/ + "mode_ac": "land", /*3 modes to determine which atmospheric correction is applied at which surfaces (default: 'land'): + 'land': SICOR (developed for land surfaces) is applied to land AND water surfaces + 'water': POLYMER (developed for water surfaces) is applied to water only + (land surfaces are not included in the L2A product) + 'combined': SICOR is applied to land and POLYMER is applied to water surfaces; + NOTE that this may result in edge effects, e.g., at coastlines*/ + "auto_download_ecmwf": false, /*FIXME this might be not needed anymore in future*/ + "enable_ice_retrieval": true, /*enable ice retrieval (default); increases accuracy of water vapour retrieval*/ "enable_cloud_screening": false, /*FIXME this is currently not implemented*/ - "scale_factor_boa_ref": 10000 /*scale factor to be applied to BOA reflectance result*/ + "scale_factor_boa_ref": 10000, /*scale factor to be applied to BOA reflectance result*/ + "threads": -1, /*number of threads for multiprocessing: + 'threads = 0': for single thread + 'threads < 0': for as many threads as there are CPUs + 'threads > 0': gives the number of threads*/ + "blocksize":100 /*block size for multiprocessing*/ + + }, "smile": { diff --git a/enpt/options/options_schema.py b/enpt/options/options_schema.py index 7a949b7f5f761110c2abc2aca40b9a5230468d76..eb9e6ee77a70c2260015f7ea04b325ebc4dc3971 100644 --- a/enpt/options/options_schema.py +++ b/enpt/options/options_schema.py @@ -82,9 +82,16 @@ enpt_schema_input = dict( atmospheric_correction=dict( type='dict', required=False, schema=dict( + polymer_root=dict(type='string', required=False), enable_ac=dict(type='boolean', required=False), + mode_ac=dict(type='string', required=False, allowed=['land', 'water', 'combined']), + auto_download_ecmwf=dict(type='boolean', required=False), + enable_ice_retrieval=dict(type='boolean', required=False), enable_cloud_screening=dict(type='boolean', required=False), scale_factor_boa_ref=dict(type='integer', required=False, min=1), + threads=dict(type='integer', required=False), + blocksize=dict(type='integer', required=False), + )), smile=dict( @@ -151,9 +158,15 @@ parameter_mapping = dict( path_reference_image=('processors', 'geometry', 'path_reference_image'), # processors > atmospheric_correction + polymer_root=('processors', 'atmospheric_correction', 'polymer_root'), enable_ac=('processors', 'atmospheric_correction', 'enable_ac'), + mode_ac=('processors', 'atmospheric_correction', 'mode_ac'), + auto_download_ecmwf=('processors', 'atmospheric_correction', 'auto_download_ecmwf'), + enable_ice_retrieval=('processors', 'atmospheric_correction', 'enable_ice_retrieval'), enable_cloud_screening=('processors', 'atmospheric_correction', 'enable_cloud_screening'), scale_factor_boa_ref=('processors', 'atmospheric_correction', 'scale_factor_boa_ref'), + threads=('processors', 'atmospheric_correction', 'threads'), + blocksize=('processors', 'atmospheric_correction', 'blocksize'), # processors > smile run_smile_P=('processors', 'smile', 'run_processor'), diff --git a/enpt/processors/atmospheric_correction/atmospheric_correction.py b/enpt/processors/atmospheric_correction/atmospheric_correction.py index 2d07f5b5dfd363eeb8ab00d84a09748b127bf178..28bf55f9679e3d6369b38708cf39ccfad2e7d52b 100644 --- a/enpt/processors/atmospheric_correction/atmospheric_correction.py +++ b/enpt/processors/atmospheric_correction/atmospheric_correction.py @@ -27,16 +27,22 @@ # You should have received a copy of the GNU Lesser General Public License along # with this program. If not, see . -"""EnPT 'atmospheric correction module. +"""EnPT atmospheric correction module. Performs the atmospheric correction of EnMAP L1B data. """ import pprint import numpy as np from multiprocessing import cpu_count +from typing import Optional +from logging import Logger from sicor.sicor_enmap import sicor_ac_enmap from sicor.options import get_options as get_ac_options +try: + from acwater.acwater import polymer_ac_enmap +except ImportError: + polymer_ac_enmap: Optional[callable] = None from ...model.images import EnMAPL1Product_SensorGeo from ...options.config import EnPTConfig @@ -51,8 +57,15 @@ class AtmosphericCorrector(object): def __init__(self, config: EnPTConfig = None): """Create an instance of AtmosphericCorrector.""" self.cfg = config + self.is_polymer_installed = polymer_ac_enmap is not None - def get_ac_options(self, enmap_ImageL1: EnMAPL1Product_SensorGeo) -> dict: + def _get_sicor_options(self, enmap_ImageL1: EnMAPL1Product_SensorGeo, land_only=False) -> dict: + """Get a dictionary containing the SICOR options. + + :param enmap_ImageL1: EnMAP Level 1 product in sensor geometry + :param land_only: True: SICOR is applied to land only; False: SICOR is applied to all pixels + :return: dictionary of SICOR options + """ path_opts = get_path_ac_options() try: @@ -65,23 +78,31 @@ class AtmosphericCorrector(object): options["retrieval"]["cpu"] = self.cfg.CPUs or cpu_count() options["retrieval"]["disable_progressbars"] = self.cfg.disable_progress_bars - return options + # temporarily disable uncertainty measures to avoid https://gitext.gfz-potsdam.de/EnMAP/sicor/-/issues/86 + options["retrieval"]["inversion"]["full"] = False + + # set land_only mode + options["retrieval"]["land_only"] = land_only except FileNotFoundError: - raise FileNotFoundError('Could not locate options file for atmospheric correction at %s.' % path_opts) + raise FileNotFoundError(f'Could not locate options file for atmospheric correction at {path_opts}') - def run_ac(self, enmap_ImageL1: EnMAPL1Product_SensorGeo) -> EnMAPL1Product_SensorGeo: - options = self.get_ac_options(enmap_ImageL1) - enmap_ImageL1.logger.debug('AC options: \n' + pprint.pformat(options)) + enmap_ImageL1.logger.debug('SICOR AC configuration: \n' + + pprint.pformat(options)) - enmap_ImageL1.set_SWIRattr_with_transformedVNIRattr('mask_landwater') + return options - # run AC - enmap_ImageL1.logger.info("Starting atmospheric correction for VNIR and SWIR detector. " - "Source radiometric unit code is '%s'." % enmap_ImageL1.meta.vnir.unitcode) + def _run_AC__land_mode(self, + enmap_ImageL1: EnMAPL1Product_SensorGeo + ) -> (np.ndarray, np.ndarray, dict): + """Run atmospheric correction in 'land' mode, i.e., use SICOR for all surfaces.""" + if 2 in enmap_ImageL1.vnir.mask_landwater[:]: + enmap_ImageL1.logger.info( + "Running atmospheric correction in 'land' mode, i.e., SICOR is applied to ALL surfaces. " + "Uncertainty is expected for water surfaces because SICOR is designed for land only.") # run SICOR - # NOTE: - enmap_l2a_vnir, enmap_l2a_swir: reflectance between 0 and 1 + # NOTE: - boa_ref_vnir, boa_ref_swir: reflectance between 0 and 1 # - res: a dictionary containing retrieval maps with path lengths of the three water phases # and several retrieval uncertainty measures # -> cwv_model, liq_model, ice_model, intercept_model, slope_model, toa_model, @@ -91,27 +112,150 @@ class AtmosphericCorrector(object): # jacobian, convergence, iterations, gain, averaging_kernel, cost_function, # dof (degrees of freedom), information_content, retrieval_noise, smoothing_error # -> SWIR geometry (?) # FIXME - enmap_l2a_vnir, enmap_l2a_swir, res = \ - sicor_ac_enmap(enmap_l1b=enmap_ImageL1, options=options, unknowns=True, logger=enmap_ImageL1.logger) + boa_ref_vnir, boa_ref_swir, land_additional_results = \ + sicor_ac_enmap(enmap_l1b=enmap_ImageL1, + options=self._get_sicor_options(enmap_ImageL1, land_only=False), + unknowns=True, + logger=enmap_ImageL1.logger) + + return boa_ref_vnir, boa_ref_swir, land_additional_results + + def _run_AC__water_mode(self, enmap_ImageL1: EnMAPL1Product_SensorGeo + ) -> (np.ndarray, np.ndarray): + """Run atmospheric correction in 'water' mode, i.e., use ACWater/Polymer for water surfaces only. + + NOTE: + - Land surfaces are NOT included in the EnMAP L2A product. + - The output radiometric unit for water surfaces is 'water leaving reflectance'. + """ + if 1 in enmap_ImageL1.vnir.mask_landwater[:]: + enmap_ImageL1.logger.info( + "Running atmospheric correction in 'water' mode, i.e., ACWater/Polymer is applied to water " + "surfaces only. Note that land surfaces will NOT be included in the EnMAP L2A product.") + + # run ACWater/Polymer for water surfaces only + # NOTE: polymer_ac_enmap() returns masked (nan) values for land + wl_ref_vnir, wl_ref_swir = \ + polymer_ac_enmap(enmap_l1b=enmap_ImageL1, + config=self.cfg, + detector='merge') + + return wl_ref_vnir, wl_ref_swir + + def _run_AC__combined_mode(self, + enmap_ImageL1: EnMAPL1Product_SensorGeo + ) -> (np.ndarray, np.ndarray, dict): + """Run atmospheric corr. in 'combined' mode, i.e., use SICOR for land and ACWater/Polymer for water surfaces. + + NOTE: + - The output radiometric units are: + - 'surface reflectance' for land surfaces + - 'water leaving reflectance' for water surfaces + - There might be visible edge effects, e.g., at coastlines. + """ + only = 'water' if 1 not in enmap_ImageL1.vnir.mask_landwater[:] else 'land' + if 1 not in enmap_ImageL1.vnir.mask_landwater[:] or \ + 2 not in enmap_ImageL1.vnir.mask_landwater[:]: + enmap_ImageL1.logger.info( + f"Running atmospheric correction in 'combined' mode, i.e., SICOR is applied to land and " + f"ACWater/Polymer is applied to water surfaces. But the input image only contains {only} surfaces.") + + # run SICOR for land surfaces only + boa_ref_vnir_land, boa_ref_swir_land, land_additional_results = \ + sicor_ac_enmap(enmap_l1b=enmap_ImageL1, + options=self._get_sicor_options(enmap_ImageL1, land_only=True), + unknowns=True, + logger=enmap_ImageL1.logger) + + # run ACWater/Polymer for water surfaces only + # NOTE: polymer_ac_enmap() returns masked (nan) values for land + wl_ref_vnir_water, wl_ref_swir_water = \ + polymer_ac_enmap(enmap_l1b=enmap_ImageL1, + config=self.cfg, + detector='merge') + + # use mask value 2 for replacing water corrected pixels + wlboa_ref_vnir = np.where((enmap_ImageL1.vnir.mask_landwater[:] == 2)[:, :, None], + wl_ref_vnir_water, + boa_ref_vnir_land) + wlboa_ref_swir = np.where((enmap_ImageL1.swir.mask_landwater[:] == 2)[:, :, None], + wl_ref_swir_water, + boa_ref_swir_land) + + return wlboa_ref_vnir, wlboa_ref_swir, land_additional_results + + @staticmethod + def _validate_AC_results(reflectance_vnir: np.ndarray, + reflectance_swir: np.ndarray, + logger: Logger): + for detectordata, detectorname in zip([reflectance_vnir, reflectance_swir], + ['VNIR', 'SWIR']): + mean0 = np.nanmean(detectordata[:, :, 0]) + std0 = np.nanstd(detectordata[:, :, 0]) + + if np.isnan(mean0) or \ + mean0 == 0 or \ + std0 == 0: + logger.warning(f'The atmospheric correction returned empty {detectorname} bands!') + + def run_ac(self, + enmap_ImageL1: EnMAPL1Product_SensorGeo + ) -> EnMAPL1Product_SensorGeo: + """Run atmospheric correction according to the specified 'mode_ac' parameter. + + :param enmap_ImageL1: input EnMAP image containing TOA reflectance (an instance EnMAPL1Product_SensorGeo) + :return: atmospherically corrected output EnMAP image containing BOA reflectance / water leaving reflectance + (an instance EnMAPL1Product_SensorGeo) + """ + enmap_ImageL1.set_SWIRattr_with_transformedVNIRattr('mask_landwater') + + enmap_ImageL1.logger.info( + f"Starting atmospheric correction for VNIR and SWIR detector in '{self.cfg.mode_ac}' mode. " + f"Source radiometric unit code is '{enmap_ImageL1.meta.vnir.unitcode}'.") + + # run the AC + if not self.is_polymer_installed and self.cfg.mode_ac in ['water', 'combined']: + # use SICOR as fallback AC for water surfaces if ACWater/Polymer is not installed + enmap_ImageL1.logger.warning( + f"The atmospheric correction mode was set to '{self.cfg.mode_ac}' but the packages ACWater/Polymer are " + f"missing. SICOR has to be used as fallback algorithm for water surfaces.") + + reflectance_vnir, reflectance_swir, land_additional_results = \ + self._run_AC__land_mode(enmap_ImageL1) + + else: + if self.cfg.mode_ac == 'combined': + reflectance_vnir, reflectance_swir, land_additional_results = \ + self._run_AC__combined_mode(enmap_ImageL1) + + elif self.cfg.mode_ac == 'water': + reflectance_vnir, reflectance_swir = \ + self._run_AC__water_mode(enmap_ImageL1) + + elif self.cfg.mode_ac == 'land': + reflectance_vnir, reflectance_swir, land_additional_results = \ + self._run_AC__land_mode(enmap_ImageL1) + + else: + raise ValueError(self.cfg.mode_ac, + "Unexpected 'mode_ac' parameter. " + "Choose one out of 'land', 'water', 'combined'.") - # validate results - for detectordata, detectorname in zip([enmap_l2a_vnir, enmap_l2a_swir], ['VNIR', 'SWIR']): - mean0, std0 = np.nanmean(detectordata[:, :, 0]), np.nanstd(detectordata[:, :, 0]) - if np.isnan(mean0) or mean0 == 0 or std0 == 0: - enmap_ImageL1.logger.warning('The atmospheric correction returned empty %s bands!' % detectorname) + # validate outputs + self._validate_AC_results(reflectance_vnir, reflectance_swir, logger=enmap_ImageL1.logger) # join results enmap_ImageL1.logger.info('Joining results of atmospheric correction.') for in_detector, out_detector in zip([enmap_ImageL1.vnir, enmap_ImageL1.swir], - [enmap_l2a_vnir, enmap_l2a_swir]): + [reflectance_vnir, reflectance_swir]): in_detector.data = (out_detector * self.cfg.scale_factor_boa_ref).astype(np.int16) # NOTE: geotransform and projection are missing due to sensor geometry in_detector.detector_meta.unit = '0-%d' % self.cfg.scale_factor_boa_ref in_detector.detector_meta.unitcode = 'BOARef' - # FIXME what about mask_clouds, mask_clouds_confidence, ac_errors? - # FIXME use cwv_model, cwc_model, toa_model also for EnPT? + # FIXME: Consider to also join SICOR's land_additional_results + # (contains three phases of water maps and several retrieval uncertainty measures) return enmap_ImageL1 diff --git a/requirements.txt b/requirements.txt index 8ff3560d800b96f5d3452c20a3f6bab322f45940..298c95f4466c8c3888a1b768ed00e87fc3d07a96 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ +git+https://gitlab.awi.de/phytooptics/acwater.git arosics>=1.0.0 cerberus geoarray>=0.9.0 diff --git a/setup.py b/setup.py index 0bf7c68986657098aed1243c263ff00b00bae450..89d746ae64a5d60cbe4d7bb6da0bf04a9edacd43 100644 --- a/setup.py +++ b/setup.py @@ -73,8 +73,9 @@ req_lint = ['flake8', 'pycodestyle', 'pydocstyle'] req_dev = req_setup + req_test + req_doc + req_lint setup( - author="Karl Segl, Daniel Scheffler, Niklas Bohn, Stéphane Guillaso", - author_email="segl@gfz-potsdam.de, danschef@gfz-potsdam.de, nbohn@gfz-potsdam.de, stephane.guillaso@gfz-potsdam.de", + author="Karl Segl, Daniel Scheffler, Niklas Bohn, Stéphane Guillaso, Brenner Silva", + author_email="segl@gfz-potsdam.de, danschef@gfz-potsdam.de, nbohn@gfz-potsdam.de, " + "stephane.guillaso@gfz-potsdam.de, brenner.silva@awi.de", classifiers=[ 'Development Status :: 4 - Beta', 'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)', diff --git a/tests/gitlab_CI_docker/context/enpt_ci.docker b/tests/gitlab_CI_docker/context/enpt_ci.docker index dc38789fc14a3796a0853fe4ca320c6b73d8a43b..a48e1434985698ee76eca868909b71e904cd22b2 100644 --- a/tests/gitlab_CI_docker/context/enpt_ci.docker +++ b/tests/gitlab_CI_docker/context/enpt_ci.docker @@ -1,14 +1,28 @@ FROM ci_base_centos:0.1 -RUN yum update -y && \ - yum install -y texlive - # copy some needed stuff to /root COPY *.yml /root/ +COPY polymer-v4.13.tar.gz /root/ # update the ci_env environment (that already contains all packages installed via 'docker_pyenvs' repo) RUN /bin/bash -i -c "\ source /root/miniconda3/bin/activate ; \ - conda update -n base -c conda-forge conda;\ + conda install -c conda-forge mamba ; \ + mamba update -n base -c conda-forge conda;\ conda activate ci_env; \ - conda env update -n ci_env -f /root/environment_enpt.yml" + mamba env update -n ci_env -f /root/environment_enpt.yml" + +# install Polymer into the ci_env environment +RUN /bin/bash -i -c "\ + source /root/miniconda3/bin/activate ci_env; \ + cd /root/ ; \ + tar -zxvf polymer-*.tar.gz ; \ + cd polymer-v4.13 ;\ + make ;\ + pip install -e . ;\ + make auxdata_common ;\ + mkdir -p ANCILLARY/ERA5/ ANCILLARY/METEO/ ;\ + mkdir -p ANCILLARY/ERA5/2017/02/18" + +# copy pre-downloaded AUX data to Polymer root to avoid at-runtime downloads +COPY era5_20170218*.nc /root/polymer-v4.13/ANCILLARY/ERA5/2017/02/18/ \ No newline at end of file diff --git a/tests/gitlab_CI_docker/context/environment_enpt.yml b/tests/gitlab_CI_docker/context/environment_enpt.yml index 6e089a50a9c6196c24db6c39f07ab1cc8e57643c..baddad52e50c9ef97f81e410ff5ebbf3bc2ff703 100644 --- a/tests/gitlab_CI_docker/context/environment_enpt.yml +++ b/tests/gitlab_CI_docker/context/environment_enpt.yml @@ -18,6 +18,15 @@ dependencies: - sensormapgeo>=0.4.0 - sicor>=0.16.0 + # Polymer AC additional conda requirements + - cdsapi + - cython + - gdal + - netcdf4 + - pygrib + - pyhdf + - xarray + - pip: - cerberus - mvgavg @@ -28,6 +37,10 @@ dependencies: - tqdm - utm + # Polymer AC additional pip requirements + - git+https://gitlab.awi.de/phytooptics/acwater.git + - ecmwf-api-client + # test/doc/lint requirements - coverage - flake8 diff --git a/tests/gitlab_CI_docker/context/environment_enpt_polymer.yml b/tests/gitlab_CI_docker/context/environment_enpt_polymer.yml new file mode 100644 index 0000000000000000000000000000000000000000..315bd9106a45d445cc4fb8cc6fc0e37e0453b99a --- /dev/null +++ b/tests/gitlab_CI_docker/context/environment_enpt_polymer.yml @@ -0,0 +1,24 @@ +name: enpt_polymer + +channels: &id1 + - conda-forge + +dependencies: + - python=3.*.* + - pip # avoids that conda uses the wrong pip + - enpt + + - cdsapi + - cython + - gdal + - netcdf4 + - numpy + - pandas + - pygrib + - pyhdf + - scipy + - xarray + + - pip: + - git+https://gitlab.awi.de/phytooptics/acwater.git + - ecmwf-api-client diff --git a/tests/test_controller.py b/tests/test_controller.py index 353cfcc292e00538f32353c9da5ad2149776d402..e5c7a4035faf09ab82881d9c44a2dfa1cd1b76b4 100644 --- a/tests/test_controller.py +++ b/tests/test_controller.py @@ -36,10 +36,11 @@ Tests for `execution.controller` module. """ from unittest import TestCase +from unittest.mock import patch import shutil from enpt.execution.controller import EnPT_Controller -from enpt.options.config import config_for_testing, config_for_testing_dlr +from enpt.options.config import config_for_testing, config_for_testing_dlr, config_for_testing_water __author__ = 'Daniel Scheffler' @@ -68,6 +69,27 @@ class Test_EnPT_Controller_DLR_testdata(TestCase): self.CTR.run_all_processors() +class Test_EnPT_Controller_DLR_testdata_ACWater(TestCase): + def setUp(self): + self.CTR = EnPT_Controller(**config_for_testing_water) + + def tearDown(self): + # NOTE: ignore_errors deletes the folder, regardless of whether it contains read-only files + shutil.rmtree(self.CTR.cfg.output_dir, ignore_errors=True) + + def test_run_all_processors(self): + self.CTR.run_all_processors() + + @patch('enpt.processors.atmospheric_correction.atmospheric_correction.polymer_ac_enmap', None) + def test_run_all_processors_without_acwater_installed(self): + """Test to run all processors while replacing polymer_ac_enmap with None using mock.patch.""" + self.CTR.run_all_processors() + + self.assertTrue("packages ACWater/Polymer are missing. " + "SICOR has to be used as fallback algorithm for water surfaces." + in self.CTR.L1_obj.logger.captured_stream) + + if __name__ == '__main__': import nose2 nose2.main() diff --git a/tests/test_l1b_reader.py b/tests/test_l1b_reader.py index fd2b38597059c9f7befdd5cf5463a9dfabb24358..d7158a90a7c29ff1718d66fe3e09d02151846904 100644 --- a/tests/test_l1b_reader.py +++ b/tests/test_l1b_reader.py @@ -189,7 +189,7 @@ class Test_L1B_Reader_DLR(unittest.TestCase): RD = L1B_Reader(config=cfg) L1_obj = RD.read_inputdata(self.tmpdir, compute_snr=False) - self.assertEquals(L1_obj.swir.detector_meta.nwvl, 130) + self.assertEqual(L1_obj.swir.detector_meta.nwvl, 130) if __name__ == '__main__':