Merge branch 'feature/separate_ac_land_water' into 'master'

Feature/separate ac land water

Closes #47

See merge request !54

AUTHORS.rst

@@ -13,6 +13,7 @@ Contributors
 * Daniel Scheffler <danschef@gfz-potsdam.de>
   (main developer of the EnPT source code)
 * Niklas Bohn <nbohn@gfz-potsdam.de>
-  (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

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)
 -------------------

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

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

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
 ***********************

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

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,

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.

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')

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": {

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'),

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 <http://www.gnu.org/licenses/>.
 
-"""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

requirements.txt

+git+https://gitlab.awi.de/phytooptics/acwater.git
 arosics>=1.0.0
 cerberus
 geoarray>=0.9.0

setup.py

@@ -73,8 +73,9 @@ req_lint = ['flake8', 'pycodestyle', 'pydocstyle']
 req_dev = req_setup + req_test + req_doc + req_lint