Added first implementation able to compute tie points between an EnMAP image...

Added first implementation able to compute tie points between an EnMAP image and a user-provided reference image.

Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>

enpt/options/config.py

@@ -130,6 +130,8 @@ config_for_testing_dlr = dict(
     # target_projection_type='Geographic',
     # target_epsg=32632,
     # target_coord_grid=[-1.37950, -1.37923, 44.60710, 44.60737],
+    enable_absolute_coreg=True,
+    path_reference_image=os.path.join(path_enptlib, '..', 'tests', 'data', 'T30TXQ_20170218T110111_B05.jp2'),
     enable_ac=True,
     enable_ice_retrieval=False,
     CPUs=32,
@@ -209,6 +211,9 @@ class EnPTConfig(object):
         :key enable_vnir_swir_coreg:
             Enable VNIR/SWIR co-registration
 
+        :key enable_absolute_coreg:
+            Enable the co-registration of the EnMAP image to the reference image given with 'path_reference_image'
+
         :key path_reference_image:
             Reference image for co-registration.
 
@@ -318,6 +323,7 @@ class EnPTConfig(object):
         # geometry
         self.enable_keystone_correction = gp('enable_keystone_correction')
         self.enable_vnir_swir_coreg = gp('enable_vnir_swir_coreg')
+        self.enable_absolute_coreg = gp('enable_absolute_coreg')
         self.path_reference_image = gp('path_reference_image')
 
         # atmospheric_correction

enpt/options/options_default.json

@@ -44,7 +44,11 @@
         "geometry": {
             "enable_keystone_correction": false,
             "enable_vnir_swir_coreg": false,
-            "path_reference_image": ""
+            "enable_absolute_coreg": false,  /*enable the co-registration of the EnMAP image to the reference image given with 'path_reference_image'*/
+            "path_reference_image": ""  /*image to be used as spatial reference for co-registration
+                                          - provide in a GDAL-readable data format
+                                          - acquisition date should be close to the one of the EnMAP scene
+                                          - cloud coverage should be low*/
         },
 
         "atmospheric_correction": {

enpt/options/options_schema.py

@@ -75,6 +75,7 @@ enpt_schema_input = dict(
                 schema=dict(
                     enable_keystone_correction=dict(type='boolean', required=False),
                     enable_vnir_swir_coreg=dict(type='boolean', required=False),
+                    enable_absolute_coreg=dict(type='boolean', required=False),
                     path_reference_image=dict(type='string', required=False),
                 )),
 
@@ -148,6 +149,7 @@ parameter_mapping = dict(
     # processors > geometry
     enable_keystone_correction=('processors', 'geometry', 'enable_keystone_correction'),
     enable_vnir_swir_coreg=('processors', 'geometry', 'enable_vnir_swir_coreg'),
+    enable_absolute_coreg=('processors', 'geometry', 'enable_absolute_coreg'),
     path_reference_image=('processors', 'geometry', 'path_reference_image'),
 
     # processors > atmospheric_correction

enpt/processors/orthorectification/orthorectification.py

@@ -58,7 +58,7 @@ class Orthorectifier(object):
         self.cfg = config
 
     @staticmethod
-    def validate_input(enmap_ImageL1):
+    def validate_input(enmap_ImageL1: EnMAPL1Product_SensorGeo):
         # check type
         if not isinstance(enmap_ImageL1, EnMAPL1Product_SensorGeo):
             raise TypeError(enmap_ImageL1, "The Orthorectifier expects an instance of 'EnMAPL1Product_SensorGeo'."

enpt/processors/spatial_optimization/__init__.py

@@ -29,4 +29,10 @@
 
 """EnPT 'spatial optimization' module."""
 
+from .spatial_optimization import Spatial_Optimizer
+
 __author__ = 'Daniel Scheffler'
+__all__ = [
+    '__author__',
+    'Spatial_Optimizer'
+]
\ No newline at end of file

enpt/processors/spatial_optimization/spatial_optimization.py

@@ -34,3 +34,70 @@ Fits the VNIR detector data to the reference image. Corrects for keystone.
 """
 
 __author__ = 'Daniel Scheffler'
+
+from arosics import COREG_LOCAL
+from geoarray import GeoArray
+
+from ...options.config import EnPTConfig
+from ...model.images.images_sensorgeo import EnMAPL1Product_SensorGeo
+from ..spatial_transform import Geometry_Transformer
+
+
+class Spatial_Optimizer(object):
+    def __init__(self, config: EnPTConfig):
+        """Create an instance of Spatial_Optimizer"""
+        self.cfg = config
+
+    def _get_enmap_band_for_matching(self,
+                                     enmap_ImageL1: EnMAPL1Product_SensorGeo)\
+            -> GeoArray:
+        """
+
+        :param enmap_ImageL1:
+        :return:
+        """
+        enmap_ImageL1.logger.warning('Statically using band 40 for co-registration.')
+        bandidx = 39
+        enmap_band_sensorgeo = enmap_ImageL1.vnir.data[:, :, bandidx]
+
+        # transform from sensor to map geometry to make ithe band usable for tie point detection
+        enmap_ImageL1.logger.info('Temporarily transforming EnMAP band %d to map geometry for co-registration.'
+                                  % (bandidx + 1))
+        GT = Geometry_Transformer(lons=enmap_ImageL1.meta.vnir.lons[:, :, bandidx],
+                                  lats=enmap_ImageL1.meta.vnir.lats[:, :, bandidx],
+                                  fill_value=0,
+                                  resamp_alg='gauss',
+                                  radius_of_influence=30,
+                                  nprocs=self.cfg.CPUs)
+
+        ref_gA = GeoArray(self.cfg.path_reference_image)
+        ref_ULx, ref_ULy = ref_gA.box.boxMapXY[0]
+
+        enmap_band_mapgeo = \
+            GeoArray(*GT.to_map_geometry(enmap_band_sensorgeo,
+                                         tgt_prj=ref_gA.prj,  # TODO correct?
+                                         tgt_coordgrid=((ref_ULx, ref_ULx + ref_gA.xgsd),
+                                                        (ref_ULy, ref_ULy + ref_gA.ygsd))
+                                         ),
+                     nodata=0)
+
+        return enmap_band_mapgeo
+
+    def _compute_tie_points(self,
+                            enmap_ImageL1: EnMAPL1Product_SensorGeo):
+        enmap_band_mapgeo = self._get_enmap_band_for_matching(enmap_ImageL1)
+
+        CRL = COREG_LOCAL(self.cfg.path_reference_image,
+                          enmap_band_mapgeo,
+                          grid_res=50)
+        TPG = CRL.tiepoint_grid
+        CRL.view_CoRegPoints(shapes2plot='vectors', hide_filtered=False)
+
+        # filter out tie points over water
+        pass  # TODO
+        a = 1  # FIXME
+
+    def optimize_geolayer(self,
+                          enmap_ImageL1: EnMAPL1Product_SensorGeo):
+        if self.cfg.enable_absolute_coreg:
+            self._compute_tie_points(enmap_ImageL1)

tests/test_spatial_optimization.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# EnPT, EnMAP Processing Tool - A Python package for pre-processing of EnMAP Level-1B data
+#
+# Copyright (C) 2019  Karl Segl (GFZ Potsdam, segl@gfz-potsdam.de), Daniel Scheffler
+# (GFZ Potsdam, danschef@gfz-potsdam.de), Niklas Bohn (GFZ Potsdam, nbohn@gfz-potsdam.de),
+# Stéphane Guillaso (GFZ Potsdam, stephane.guillaso@gfz-potsdam.de)
+#
+# This software was developed within the context of the EnMAP project supported
+# by the DLR Space Administration with funds of the German Federal Ministry of
+# Economic Affairs and Energy (on the basis of a decision by the German Bundestag:
+# 50 EE 1529) and contributions from DLR, GFZ and OHB System AG.
+#
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version. Please note the following exception: `EnPT` depends on tqdm, which
+# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files
+# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore".
+# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+# details.
+#
+# 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/>.
+
+"""
+test_spatial_optimization
+-------------------------
+
+Tests for `processors.spatial_optimization.spatial_optimization` module.
+"""
+
+import os
+from unittest import TestCase
+from zipfile import ZipFile
+import tempfile
+import shutil
+
+from enpt.processors.spatial_optimization import Spatial_Optimizer
+from enpt.options.config import config_for_testing, config_for_testing_dlr, EnPTConfig
+from enpt.io.reader import L1B_Reader
+
+__author__ = 'Daniel Scheffler'
+
+
+class Test_Spatial_Optimizer(TestCase):
+    def setUp(self):
+        self.config = EnPTConfig(**config_for_testing)
+
+        # create a temporary directory
+        # NOTE: This must exist during the whole runtime of Test_Spatial_Optimizer, otherwise
+        #       Spatial_Optimizer.optimize_geolayer will fail to read some files.
+        self.tmpdir = tempfile.mkdtemp(dir=self.config.working_dir)
+
+        # get EnMAP L1 object in sensor geometry
+        with ZipFile(self.config.path_l1b_enmap_image, "r") as zf:
+            zf.extractall(self.tmpdir)
+            self.L1_obj = L1B_Reader(config=self.config).read_inputdata(
+                root_dir_main=os.path.join(self.tmpdir,
+                                           os.path.splitext(os.path.basename(self.config.path_l1b_enmap_image))[0]),
+                compute_snr=False)
+
+    def tearDown(self):
+        shutil.rmtree(self.tmpdir)
+
+    def test_optimize_geolayer(self):
+        SO = Spatial_Optimizer(config=self.config)
+        out = SO.optimize_geolayer(self.L1_obj)
+
+        # self.assertIsInstance(L2_obj, EnMAPL2Product_MapGeo)
+        # self.assertTrue(L2_obj.data.is_map_geo)
+        # self.assertGreater(L2_obj.data.shape[0], self.L1_obj.vnir.data.shape[0])
+        # self.assertNotEqual(L2_obj.data.shape[1], self.L1_obj.vnir.data.shape[1])
+        # self.assertEqual(L2_obj.data.ndim, self.L1_obj.vnir.data.ndim)
+        # self.assertTrue(np.isclose(np.mean(self.L1_obj.vnir.data[:, :, 0]),
+        #                            np.mean(L2_obj.data[:, :, 0][L2_obj.data[:, :, 0] != L2_obj.data.nodata]),
+        #                            rtol=0.01
+        #                            ))
+
+
+class Test_Orthorectifier_DLR(TestCase):
+    def setUp(self):
+        self.config = EnPTConfig(**config_for_testing_dlr)
+
+        # create a temporary directory
+        # NOTE: This must exist during the whole runtime of Test_Spatial_Optimizer, otherwise
+        #       Spatial_Optimizer.optimize_geolayer will fail to read some files.
+        self.tmpdir = tempfile.mkdtemp(dir=self.config.working_dir)
+
+        # get EnMAP L1 object in sensor geometry
+        with ZipFile(self.config.path_l1b_enmap_image, "r") as zf:
+            zf.extractall(self.tmpdir)
+            self.L1_obj = L1B_Reader(config=self.config).read_inputdata(
+                root_dir_main=self.tmpdir,
+                compute_snr=False)
+
+    def tearDown(self):
+        shutil.rmtree(self.tmpdir)
+
+    def test_optimize_geolayer(self):
+        SO = Spatial_Optimizer(config=self.config)
+        out = SO.optimize_geolayer(self.L1_obj)
+
+        # self.assertIsInstance(L2_obj, EnMAPL2Product_MapGeo)
+        # self.assertTrue(L2_obj.data.is_map_geo)
+        # self.assertGreater(L2_obj.data.shape[0], self.L1_obj.vnir.data.shape[0])
+        # self.assertNotEqual(L2_obj.data.shape[1], self.L1_obj.vnir.data.shape[1])
+        # self.assertEqual(L2_obj.data.ndim, self.L1_obj.vnir.data.ndim)
+        # self.assertTrue(np.isclose(np.mean(self.L1_obj.vnir.data[:, :, 0]),
+        #                            np.mean(L2_obj.data[:, :, 0][L2_obj.data[:, :, 0] != L2_obj.data.nodata]),
+        #                            rtol=0.01
+        #                            ))