Commit 8542621a authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Added basic compatibility with images that have a rotation in the map info...

Added basic compatibility with images that have a rotation in the map info (fixes #60

). Fixed incorrect tolerance in COREG.equalize_pixGrids() which speeds up COREG_LOCAL a lot. Bumped version.
Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent 41b2c1f2
Pipeline #25531 passed with stages
in 5 minutes and 7 seconds
......@@ -2,11 +2,13 @@
History
=======
1.4.9 (coming soon)
1.5.0 (2021-07-13)
-------------------
* COREG.show_matchWin() now also works when no valid shift was found yet.
* Updated minimal version of geoarray to fix a sporadic TypeError when writing the coregistered result.
* Added basic compatibility with images that have a rotation in the map info (fixes #60).
* Fixed incorrect tolerance in COREG.equalize_pixGrids() which speeds up COREG_LOCAL a lot.
1.4.8 (2021-07-06)
......
......@@ -407,6 +407,7 @@ class COREG(object):
# try:
gdal.AllRegister()
self._check_and_handle_metaRotation()
self._get_image_params()
self._set_outpathes(im_ref, im_tgt)
self.grid2use = 'ref' if self.shift.xgsd <= self.ref.xgsd else 'shift'
......@@ -576,26 +577,55 @@ class COREG(object):
assert px_covered > 16 * 16, \
'Overlap area covers only %s pixels. At least 16*16 pixels are needed.' % px_covered
def _check_and_handle_metaRotation(self):
"""Check if the provided input data have a metadata rotation and if yes, correct it.
In case there is a rotation, the GDAL GeoTransform is not 0 at positions 2 or 4. So far, AROSICS does not
handle such rotations, so the resampling is needed to make things work.
"""
gA_ref = GeoArray(self.params['im_ref'])
gA_tgt = GeoArray(self.params['im_tgt'])
msg = 'The %s image needs to be resampled because it has a row/column rotation in ' \
'its map info which is not handled by AROSICS.'
if GEO.has_metaRotation(gA_ref):
warnings.warn(msg % 'reference')
self.params['im_ref'] = GEO.remove_metaRotation(gA_ref)
if GEO.has_metaRotation(gA_tgt):
warnings.warn(msg % 'target')
self.params['im_tgt'] = GEO.remove_metaRotation(gA_tgt)
@property
def are_pixGrids_equal(self):
return prj_equal(self.ref.prj, self.shift.prj) and \
is_coord_grid_equal(self.ref.gt, *self.shift.xygrid_specs, tolerance=1e-8)
def equalize_pixGrids(self) -> None:
"""Equalize image grids and projections of reference and target image (align target to reference)."""
if not (prj_equal(self.ref.prj, self.shift.prj) and
is_coord_grid_equal(self.ref.gt, *self.shift.xygrid_specs, tolerance=1e8)):
"""Equalize image grids and projections of reference and target image (align target to reference).
NOTE: This method is only called by COREG_LOCAL to speed up processing during detection of displacements.
"""
if not self.are_pixGrids_equal or GEO.has_metaRotation(self.ref) or GEO.has_metaRotation(self.shift):
if not self.q:
print("Equalizing pixel grids and projections of reference and target image...")
def equalize(gA_from: GeoArray, gA_to: GeoArray) -> GeoArray:
if gA_from.bands > 1:
gA_from = gA_from.get_subset(zslice=slice(gA_from.band4match, gA_from.band4match + 1))
gA_from.reproject_to_new_grid(prototype=gA_to, CPUs=self.CPUs)
gA_from.band4match = 0 # after resampling there is only one band in the GeoArray
return gA_from
if self.grid2use == 'ref':
# resample target image to reference image
if self.shift.bands > 1:
self.shift = self.shift.get_subset(zslice=slice(self.shift.band4match, self.shift.band4match + 1))
self.shift.reproject_to_new_grid(prototype=self.ref, CPUs=self.CPUs)
self.shift.band4match = 0 # after resampling there is only one band in the GeoArray
# resample target to reference image
self.shift = equalize(gA_from=self.shift, gA_to=self.ref)
else:
# resample reference image to target image
if self.ref.bands > 1:
self.ref = self.ref.get_subset(zslice=slice(self.ref.band4match, self.ref.band4match + 1))
self.ref.reproject_to_new_grid(prototype=self.shift, CPUs=self.CPUs)
self.ref.band4match = 0 # after resampling there is only one band in the GeoArray
# resample reference to target image
self.ref = equalize(gA_from=self.ref, gA_to=self.shift)
# self.ref.gt = (self.ref.gt[0], 1, self.ref.gt[2], self.ref.gt[3], self.ref.gt[4], -1)
# self.shift.gt = (self.shift.gt[0], 1, self.shift.gt[2], self.shift.gt[3], self.shift.gt[4], -1)
......
......@@ -41,6 +41,7 @@ from geopandas import GeoDataFrame # noqa F401
from .Tie_Point_Grid import Tie_Point_Grid
from .CoReg import COREG
from .DeShifter import DESHIFTER
from .geometry import has_metaRotation, remove_metaRotation
from py_tools_ds.geo.coord_trafo import transform_any_prj, reproject_shapelyGeometry
from py_tools_ds.geo.map_info import geotransform2mapinfo
from geoarray import GeoArray
......@@ -268,6 +269,9 @@ class COREG_LOCAL(object):
self.params = dict([x for x in locals().items() if x[0] != "self" and not x[0].startswith('__')])
# NOTE: self.imref and self.im2shift are handled completely independent from self.COREG_obj.ref and
# self.COREG_obj.shift. self.COREG_obj.ref and self.COREG_obj.shift are used for shift calculation and
# correction is applied to self.im2shift.
self.imref = GeoArray(im_ref, nodata=nodata[0], progress=progress, q=q)
self.im2shift = GeoArray(im_tgt, nodata=nodata[1], progress=progress, q=q)
self.path_out = path_out # updated by self.set_outpathes
......@@ -315,6 +319,9 @@ class COREG_LOCAL(object):
gdal.AllRegister()
# resample input data in case there is a metadata rotation (not handled by AROSICS)
self._check_and_handle_metaRotation()
try:
# ignore_errors must be False because in case COREG init fails, coregistration for the whole scene fails
self.COREG_obj = COREG(self.imref, self.im2shift,
......@@ -360,6 +367,40 @@ class COREG_LOCAL(object):
self.deshift_results = None # set by self.correct_shifts()
self._success = None # set by self.success property
def _check_and_handle_metaRotation(self):
"""Check if the provided input data have a metadata rotation and if yes, correct it AND equalize grids.
In case there is a rotation, the GDAL GeoTransform is not 0 at positions 2 or 4. So far, AROSICS does not
handle such rotations, so the resampling is needed to make things work. The pixel grid equalization is also
done here to avoid a double-resampling (grid would be equalized by COREG.equalize_pixGrids() otherwise).
"""
grid2use = 'ref' if self.im2shift.xgsd <= self.imref.xgsd else 'shift'
if has_metaRotation(self.imref) or has_metaRotation(self.im2shift):
msg = 'The %s image needs to be resampled because it has a row/column rotation in '\
'its map info which is not handled by AROSICS.'
if grid2use == 'ref':
if has_metaRotation(self.imref):
warnings.warn(msg % 'reference')
self.imref = remove_metaRotation(self.imref)
# resample target to reference image
if not self.q:
print('Adapting the target image pixel grid to the one of the reference image for shift detection.')
self.im2shift.reproject_to_new_grid(prototype=self.imref, CPUs=self.CPUs)
else:
# remove any metadata rotation (a rotation that only exists in the map info)
if has_metaRotation(self.im2shift):
warnings.warn(msg % 'target')
self.im2shift = remove_metaRotation(self.im2shift)
# resample reference to target image
print('Adapting the reference image pixel grid to the one of the target image for shift detection.')
self.imref.reproject_to_new_grid(prototype=self.im2shift, CPUs=self.CPUs)
def check_if_fftw_works(self) -> None:
"""Assign the attribute 'fftw_works' to self.COREG_obj by executing shift calculation once with muted output."""
# calculate global shift once in order to check is fftw works
......@@ -757,7 +798,16 @@ class COREG_LOCAL(object):
if max_GCP_count:
self.coreg_info['GCPList'] = self.coreg_info['GCPList'][:max_GCP_count]
DS = DESHIFTER(self.im2shift, self.coreg_info,
# make sure the correction is applied to the original target image
im2shift = GeoArray(self.params['im_tgt'], nodata=self.nodata[1], progress=self.progress, q=self.q)
if has_metaRotation(im2shift):
# resample the target image because (so far) the computed shifts cannot be applied to a dataset with
# a metadata rotation (GDAL GeoTransform not 0 at positons 2 and 4)
im2shift = remove_metaRotation(im2shift)
# apply the correction
DS = DESHIFTER(im2shift, self.coreg_info,
path_out=self.path_out,
fmt_out=self.fmt_out,
out_crea_options=self.out_creaOpt,
......
......@@ -360,6 +360,8 @@ class Tie_Point_Grid(object):
# NOTE: actually grid res should be also changed here because self.shift.xgsd changes and grid res is
# connected to that
self.COREG_obj.equalize_pixGrids()
self.ref = self.COREG_obj.ref
self.shift = self.COREG_obj.shift
# validate reference and target image inputs
assert self.ref.footprint_poly # this also checks for mask_nodata and nodata value
......
......@@ -23,6 +23,8 @@
import warnings
import sys
import os
from typing import Union
# custom
import numpy as np
......@@ -31,6 +33,7 @@ from geopandas import GeoDataFrame
# internal modules
from py_tools_ds.geo.coord_calc import calc_FullDataset_corner_positions
from py_tools_ds.geo.coord_trafo import pixelToMapYX, imYX2mapYX
from py_tools_ds.geo.raster.reproject import warp_ndarray
from geoarray import GeoArray
__author__ = 'Daniel Scheffler'
......@@ -133,3 +136,23 @@ def get_GeoArrayPosition_from_boxImYX(boxImYX):
rS, cS = boxImYX[0] # UL
rE, cE = boxImYX[2] # LR
return rS, rE - 1, cS, cE - 1 # -1 because boxImYX represents outer box and includes the LR corner of LR pixel
def has_metaRotation(path_or_geoarray: Union[GeoArray, str]):
"""Return True if there is a row or column rotation due to the given GDAL GeoTransform tuple."""
gt = GeoArray(path_or_geoarray).gt
return gt[2] or gt[4]
def remove_metaRotation(gA_rot: GeoArray, rspAlg='cubic') -> GeoArray:
"""Remove any metadata rotation (a rotation that only exists in the map info)."""
gA = GeoArray(*warp_ndarray(gA_rot[:], gA_rot.gt, gA_rot.prj,
rspAlg=rspAlg,
# out_gsd=(gA_rot.xgsd, gA_rot.ygsd)
),
nodata=gA_rot.nodata)
gA.basename = os.path.basename(gA.basename)
gA.meta = gA.meta
return gA
......@@ -22,5 +22,5 @@
# with this program. If not, see <http://www.gnu.org/licenses/>.
__version__ = '1.4.8'
__versionalias__ = '2021-07-06_01'
__version__ = '1.5.0'
__versionalias__ = '2021-07-13_01'
......@@ -167,6 +167,30 @@ class CompleteWorkflow_INTER1_S2A_S2A(unittest.TestCase):
footprint_poly_tgt=None))
self.assertTrue(CR.success)
def test_shift_calculation_with_metaRotation(self):
"""Test with default parameters - should compute X/Y shifts properly and write the de-shifted target image."""
# overwrite gt and prj
ref = GeoArray(self.ref_path)
ref.to_mem()
ref.filePath = None
ref.gt = [330000, 10, 0.0, 5862000, 0.0, -10]
tgt = GeoArray(self.tgt_path)
tgt.to_mem()
tgt.filePath = None
# tgt.gt = [335440, 5.8932, 0.0, 5866490, 0.0, -10.1]
tgt.gt = [335440, 10, 0.00001, 5866490, 0.00001, -10]
CR = self.run_shift_detection_correction(ref, tgt,
**dict(self.coreg_kwargs,
# ws=(512, 512),
wp=(341500.0, 5861440.0),
footprint_poly_ref=None,
footprint_poly_tgt=None,
max_shift=35))
CR.show_matchWin(interactive=False, after_correction=None)
self.assertTrue(CR.success)
def test_shift_calculation_inmem_gAs_path_out_auto(self):
"""Test input parameter path_out='auto' in case input reference/ target image are in-memory GeoArrays."""
ref = GeoArray(np.random.randint(1, 100, (1000, 1000)))
......
......@@ -175,6 +175,28 @@ class CompleteWorkflow_INTER1_S2A_S2A(unittest.TestCase):
CRL.calculate_spatial_shifts()
# CRL.view_CoRegPoints()
def test_calculation_of_tie_point_grid_with_metaRotation(self):
"""Test with default parameters - should compute X/Y shifts properly and write the de-shifted target image."""
# overwrite gt and prj
ref = GeoArray(self.ref_path)
ref.to_mem()
ref.filePath = None
ref.gt = [330000, 10, 0.00001, 5862000, 0.00001, -10]
tgt = GeoArray(self.tgt_path)
tgt.to_mem()
tgt.filePath = None
# tgt.gt = [335440, 5.8932, 0.0, 5866490, 0.0, -10.1]
tgt.gt = [335440, 10, 0.00001, 5866490, 0.00001, -10]
# get instance of COREG_LOCAL object
CRL = COREG_LOCAL(ref, tgt, **dict(CPUs=32,
**self.coreg_kwargs))
CRL.calculate_spatial_shifts()
# CRL.view_CoRegPoints()
self.assertTrue(CRL.success)
if __name__ == '__main__':
import nose2
......
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