Commit 069d0a27 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'feature/add_coregistration' into 'master'

Feature/add coregistration

Closes #33

See merge request !49
parents c6ae228e 9d00e658
Pipeline #12688 passed with stages
in 8 minutes and 27 seconds
*.zip filter=lfs diff=lfs merge=lfs -text
*.bsq filter=lfs diff=lfs merge=lfs -text
*.tif filter=lfs diff=lfs merge=lfs -text
......@@ -2,6 +2,14 @@
History
=======
0.15.0 (2020-09-21)
-------------------
* Added functionality to apply co-registration between an EnMAP image and a user-provided spatial reference dataset
(still needs to be improved but can already be used). This includes: Spatial_Optimizer class, Test_Spatial_Optimizer
class, updated config parameters, spatial reference test image.
0.14.1 (2020-09-01)
-------------------
......
......@@ -80,6 +80,8 @@ This software was developed within the context of the EnMAP project supported by
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.
Sentinel-2 spatial reference test data have been provided by ESA.
This package was created with Cookiecutter_ and the `audreyr/cookiecutter-pypackage`_ project template.
.. _Cookiecutter: https://github.com/audreyr/cookiecutter
......
......@@ -128,8 +128,13 @@ class EnPT_Controller(object):
def run_dem_processor(self):
self.L1_obj.get_preprocessed_dem()
def run_geometry_processor(self):
pass
def run_spatial_optimization(self):
# get a new instance of radiometric transformer
from ..processors.spatial_optimization import Spatial_Optimizer
SpO = Spatial_Optimizer(self.cfg)
# run optimization
self.L1_obj = SpO.optimize_geolayer(self.L1_obj)
def run_atmospheric_correction(self):
"""Run atmospheric correction only."""
......@@ -157,7 +162,9 @@ class EnPT_Controller(object):
self.read_L1B_data()
if self.cfg.run_deadpix_P:
self.L1_obj.correct_dead_pixels()
# self.run_toaRad2toaRef() # this is only needed for geometry processor but AC expects radiance
if self.cfg.enable_absolute_coreg:
# self.run_toaRad2toaRef() # this is only needed for geometry processor but AC expects radiance
self.run_spatial_optimization()
self.run_dem_processor()
if self.cfg.enable_ac:
self.run_atmospheric_correction()
......@@ -170,7 +177,7 @@ class EnPT_Controller(object):
self.L1_obj.logger.info('Skipping atmospheric correction as configured and '
'computing top-of-atmosphere reflectance instead.')
self.run_toaRad2toaRef()
self.run_geometry_processor()
# self.run_spatial_optimization()
self.run_orthorectification()
self.write_output()
......
......@@ -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__sub.tif'),
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
......
......@@ -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": {
......
......@@ -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
......
......@@ -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'."
......
......@@ -29,4 +29,10 @@
"""EnPT 'spatial optimization' module."""
from .spatial_optimization import Spatial_Optimizer
__author__ = 'Daniel Scheffler'
__all__ = [
'__author__',
'Spatial_Optimizer'
]
......@@ -34,3 +34,206 @@ Fits the VNIR detector data to the reference image. Corrects for keystone.
"""
__author__ = 'Daniel Scheffler'
import numpy as np
from typing import Optional
from arosics import COREG_LOCAL
from geoarray import GeoArray
from py_tools_ds.geo.coord_trafo import reproject_shapelyGeometry, transform_coordArray
from py_tools_ds.geo.projection import EPSG2WKT
from ...options.config import EnPTConfig
from ...model.images.images_sensorgeo import EnMAPL1Product_SensorGeo
from ..spatial_transform import Geometry_Transformer, Geometry_Transformer_3D
class Spatial_Optimizer(object):
def __init__(self, config: EnPTConfig):
"""Create an instance of Spatial_Optimizer"""
self.cfg = config
self._ref_Im: Optional[GeoArray, None] = GeoArray(self.cfg.path_reference_image)
self._EnMAP_Im: Optional[EnMAPL1Product_SensorGeo, None] = None
self._EnMAP_band: Optional[GeoArray, None] = None
self._EnMAP_mask: Optional[GeoArray, None] = None
def _get_enmap_band_for_matching(self)\
-> GeoArray:
"""Return the EnMAP band to be used in co-registration in the projection of the reference image."""
self._EnMAP_Im.logger.warning('Statically using band 40 for co-registration.')
bandidx = 39 # FIXME hardcoded
enmap_band_sensorgeo = self._EnMAP_Im.vnir.data[:, :, bandidx]
# transform from sensor to map geometry to make it usable for tie point detection
self._EnMAP_Im.logger.info('Temporarily transforming EnMAP band %d to map geometry for co-registration.'
% (bandidx + 1))
GT = Geometry_Transformer(lons=self._EnMAP_Im.meta.vnir.lons[:, :, bandidx],
lats=self._EnMAP_Im.meta.vnir.lats[:, :, bandidx],
fill_value=0,
resamp_alg='gauss',
radius_of_influence=30,
nprocs=self.cfg.CPUs)
self._EnMAP_band = \
GeoArray(*GT.to_map_geometry(enmap_band_sensorgeo,
tgt_prj=self._ref_Im.prj, # TODO correct?
tgt_coordgrid=self._ref_Im.xygrid_specs),
nodata=0)
return self._EnMAP_band
def _get_enmap_mask_for_matching(self)\
-> GeoArray:
"""Return the EnMAP mask to be used in co-registration in the projection of the reference image."""
# use the water mask
enmap_mask_sensorgeo = self._EnMAP_Im.vnir.mask_landwater[:] == 2 # 2 is water here
# transform from sensor to map geometry to make it usable for tie point detection
self._EnMAP_Im.logger.info('Temporarily transforming EnMAP water mask to map geometry for co-registration.')
GT = Geometry_Transformer(lons=self._EnMAP_Im.meta.vnir.lons[:, :, 39], # FIXME hardcoded
lats=self._EnMAP_Im.meta.vnir.lats[:, :, 39],
fill_value=0,
resamp_alg='nearest',
nprocs=self.cfg.CPUs)
self._EnMAP_mask = \
GeoArray(*GT.to_map_geometry(enmap_mask_sensorgeo,
tgt_prj=self._ref_Im.prj, # TODO correct?
tgt_coordgrid=self._ref_Im.xygrid_specs),
nodata=0)
return self._EnMAP_mask
def _compute_tie_points(self):
CRL = COREG_LOCAL(self.cfg.path_reference_image,
self._EnMAP_band,
grid_res=40,
max_shift=5,
nodata=(self._ref_Im.nodata, 0),
footprint_poly_tgt=reproject_shapelyGeometry(self._EnMAP_Im.meta.vnir.ll_mapPoly,
4326, self._EnMAP_band.epsg),
mask_baddata_tgt=self._EnMAP_mask
)
TPG = CRL.tiepoint_grid
# CRL.view_CoRegPoints(shapes2plot='vectors', hide_filtered=False, figsize=(20, 20),
# savefigPath='/home/gfz-fe/scheffler/temp/EnPT/Archachon_AROSICS_tiepoints.png')
valid_tiepoints = TPG.CoRegPoints_table[TPG.CoRegPoints_table.OUTLIER.__eq__(False)].copy()
return valid_tiepoints
@staticmethod
def _interpolate_tiepoints_into_space(tiepoints, outshape, metric='ABS_SHIFT'):
rows = np.array(tiepoints.Y_IM)
cols = np.array(tiepoints.X_IM)
data = np.array(tiepoints[metric])
from time import time
t0 = time()
# https://github.com/agile-geoscience/xlines/blob/master/notebooks/11_Gridding_map_data.ipynb
from scipy.interpolate import Rbf
# f = Rbf(cols, rows, data, function='linear')
# f = Rbf(cols, rows, data)
# data_full = f(*np.meshgrid(np.arange(outshape[1]),
# np.arange(outshape[0])))
# rows_lowres = np.arange(0, outshape[0] + 10, 10)
# cols_lowres = np.arange(0, outshape[1] + 10, 10)
rows_lowres = np.arange(0, outshape[0] + 5, 5)
cols_lowres = np.arange(0, outshape[1] + 5, 5)
f = Rbf(cols, rows, data)
data_interp_lowres = f(*np.meshgrid(cols_lowres, rows_lowres))
# https://stackoverflow.com/questions/24978052/interpolation-over-regular-grid-in-python
# from sklearn.gaussian_process import GaussianProcess
# gp = GaussianProcess(theta0=0.1, thetaL=.001, thetaU=1., nugget=0.01)
# gp.fit(X=np.column_stack([rr[vals], cc[vals]]), y=M[vals])
# rr_cc_as_cols = np.column_stack([rr.flatten(), cc.flatten()])
# interpolated = gp.predict(rr_cc_as_cols).reshape(M.shape)
from scipy.interpolate import RegularGridInterpolator
RGI = RegularGridInterpolator(points=[cols_lowres, rows_lowres],
values=data_interp_lowres.T, # must be in shape [x, y]
method='linear',
bounds_error=False)
rows_full = np.arange(outshape[0])
cols_full = np.arange(outshape[1])
data_full = RGI(np.dstack(np.meshgrid(cols_full, rows_full)))
print('interpolation runtime: %.2fs' % (time() - t0))
# from matplotlib import pyplot as plt
# plt.figure()
# im = plt.imshow(data_full)
# plt.colorbar(im)
# plt.scatter(cols, rows, c=data, edgecolors='black')
# plt.title(metric)
# plt.show()
return data_full
def optimize_geolayer(self,
enmap_ImageL1: EnMAPL1Product_SensorGeo):
self._EnMAP_Im = enmap_ImageL1
self._get_enmap_band_for_matching()
self._get_enmap_mask_for_matching()
enmap_ImageL1.logger.info('Computing tie points between the EnMAP image and the given spatial reference image.')
tiepoints = self._compute_tie_points()
enmap_ImageL1.logger.info('Generating misregistration array.')
xshift_map = self._interpolate_tiepoints_into_space(tiepoints,
self._EnMAP_band.shape,
metric='X_SHIFT_M')
yshift_map = self._interpolate_tiepoints_into_space(tiepoints,
self._EnMAP_band.shape,
metric='Y_SHIFT_M')
ULx, ULy = self._EnMAP_band.box.boxMapXY[0]
xgsd, ygsd = self._EnMAP_band.xgsd, self._EnMAP_band.ygsd
rows, cols = self._EnMAP_band.shape
xgrid_map, ygrid_map = np.meshgrid(np.arange(ULx, ULx + cols * xgsd, xgsd),
np.arange(ULy, ULy - rows * ygsd, -ygsd))
xgrid_map_coreg = xgrid_map + xshift_map
ygrid_map_coreg = ygrid_map + yshift_map
# transform map to sensor geometry
enmap_ImageL1.logger.info('Transforming spatial optimization results back to sensor geometry.')
lons_band = self._EnMAP_Im.meta.vnir.lons[:, :, 39] # FIXME hardcoded
lats_band = self._EnMAP_Im.meta.vnir.lats[:, :, 39]
GT = Geometry_Transformer_3D(lons=np.repeat(lons_band[:, :, np.newaxis], 2, axis=2),
lats=np.repeat(lats_band[:, :, np.newaxis], 2, axis=2),
fill_value=0,
# resamp_alg='bilinear',
resamp_alg='gauss',
nprocs=self.cfg.CPUs)
geolayer_sensorgeo = \
GT.to_sensor_geometry(GeoArray(np.dstack([xgrid_map_coreg,
ygrid_map_coreg]),
geotransform=self._EnMAP_band.gt,
projection=self._EnMAP_band.prj))
enmap_ImageL1.logger.info('Applying results of spatial optimization to geolayer.')
lons_coreg, lats_coreg = transform_coordArray(prj_src=self._ref_Im.prj,
prj_tgt=EPSG2WKT(4326),
Xarr=geolayer_sensorgeo[:, :, 0],
Yarr=geolayer_sensorgeo[:, :, 1])
diffs_lons_coreg = lons_band - lons_coreg
diffs_lats_coreg = lats_band - lats_coreg
# enmap_ImageL1.meta.vnir.lons -= diffs_lons_coreg[:, :, np.newaxis]
# enmap_ImageL1.meta.vnir.lats -= diffs_lats_coreg[:, :, np.newaxis]
# enmap_ImageL1.meta.swir.lons -= diffs_lons_coreg[:, :, np.newaxis]
# enmap_ImageL1.meta.swir.lats -= diffs_lats_coreg[:, :, np.newaxis]
enmap_ImageL1.meta.vnir.lons = enmap_ImageL1.meta.vnir.lons - diffs_lons_coreg[:, :, np.newaxis]
enmap_ImageL1.meta.vnir.lats = enmap_ImageL1.meta.vnir.lats - diffs_lats_coreg[:, :, np.newaxis]
enmap_ImageL1.meta.swir.lons = enmap_ImageL1.meta.swir.lons - diffs_lons_coreg[:, :, np.newaxis]
enmap_ImageL1.meta.swir.lats = enmap_ImageL1.meta.swir.lats - diffs_lats_coreg[:, :, np.newaxis]
return enmap_ImageL1
......@@ -27,6 +27,6 @@
# 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/>.
__version__ = '0.14.1'
__versionalias__ = '20200901.01'
__version__ = '0.15.0'
__versionalias__ = '20200921.01'
__author__ = 'Daniel Scheffler'
......@@ -44,6 +44,7 @@ dependencies:
- geoarray>=0.8.11
- py_tools_ds>=0.14.25
- sensormapgeo>=0.4.0
- arosics>0.9.2
- cerberus
- jsmin
- tqdm
......
#!/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.
"""
from unittest import TestCase
from zipfile import ZipFile
import tempfile
import shutil
import numpy as np
from enpt.processors.spatial_optimization import Spatial_Optimizer
from enpt.options.config import config_for_testing_dlr, EnPTConfig
from enpt.io.reader import L1B_Reader
from enpt.model.images.images_sensorgeo import EnMAPL1Product_SensorGeo
__author__ = 'Daniel Scheffler'
class Test_Spatial_Optimizer(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)
L1_obj = SO.optimize_geolayer(self.L1_obj)
self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo)
self.assertNotEqual(np.mean(L1_obj.meta.vnir.lons), 0)
self.assertNotEqual(np.std(L1_obj.meta.vnir.lons), 0)
self.assertNotEqual(np.mean(L1_obj.meta.vnir.lats), 0)
self.assertNotEqual(np.std(L1_obj.meta.vnir.lats), 0)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment