Commit f4db0513 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Fix for warping artifacts or empty output in case only very few tie points could be identified.

updated __version__ and __versionalias__
parent 95fb0774
......@@ -3,6 +3,7 @@ __author__='Daniel Scheffler'
import warnings
import os
from copy import copy
# custom
try:
......@@ -21,6 +22,7 @@ from .Tie_Point_Grid import Tie_Point_Grid
from .CoReg import COREG
from .DeShifter import DESHIFTER
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
......@@ -291,6 +293,8 @@ class COREG_LOCAL(object):
progress = self.progress,
v = self.v,
q = self.q)
self._tiepoint_grid.get_CoRegPoints_table()
if self.v:
print('Visualizing CoReg points grid...')
self.view_CoRegPoints(figsize=(10,10))
......@@ -496,6 +500,16 @@ class COREG_LOCAL(object):
return map_osm
def _get_updated_map_info_meanShifts(self):
"""Returns the updated map info of the target image, shifted on the basis of the mean X/Y shifts."""
original_map_info = geotransform2mapinfo(self.im2shift.gt, self.im2shift.prj)
updated_map_info = copy(original_map_info)
updated_map_info[3] = str(float(original_map_info[3]) + self.tiepoint_grid.mean_x_shift_map)
updated_map_info[4] = str(float(original_map_info[4]) + self.tiepoint_grid.mean_y_shift_map)
return updated_map_info
@property
def coreg_info(self):
if self._coreg_info:
......@@ -503,6 +517,12 @@ class COREG_LOCAL(object):
else:
self._coreg_info = {
'GCPList' : self.tiepoint_grid.GCPList,
'mean_shifts_px' : {'x':self.tiepoint_grid.mean_x_shift_px,
'y':self.tiepoint_grid.mean_y_shift_px},
'mean_shifts_map' : {'x':self.tiepoint_grid.mean_x_shift_map,
'y':self.tiepoint_grid.mean_y_shift_map},
'updated map info means': self._get_updated_map_info_meanShifts(),
'original map info' : geotransform2mapinfo(self.imref.gt, self.imref.prj),
'reference projection' : self.imref.prj,
'reference geotransform': self.imref.gt,
'reference grid' : [ [self.imref.gt[0], self.imref.gt[0]+self.imref.gt[1]],
......@@ -513,12 +533,15 @@ class COREG_LOCAL(object):
return self.coreg_info
def correct_shifts(self, max_GCP_count=None, cliptoextent=False):
def correct_shifts(self, max_GCP_count=None, cliptoextent=False, min_points_local_corr=5):
"""Performs a local shift correction using all points from the previously calculated geometric quality grid
that contain valid matches as GCP points.
:param max_GCP_count: <int> maximum number of GCPs to use
:param cliptoextent: <bool> whether to clip the output image to its real extent
:param min_points_local_corr: <int> number of valid tie points, below which a global shift correction is performed
instead of a local correction (global X/Y shift is then computed as the mean shift
of the remaining points)(default: 5 tie points)
:return:
"""
......@@ -529,19 +552,20 @@ class COREG_LOCAL(object):
coreg_info['GCPList'] = coreg_info['GCPList'][:max_GCP_count]
DS = DESHIFTER(self.im2shift, coreg_info,
path_out = self.path_out,
fmt_out = self.fmt_out,
out_crea_options = self.out_creaOpt,
align_grids = self.align_grids,
match_gsd = self.match_gsd,
out_gsd = self.out_gsd,
target_xyGrid = self.target_xyGrid,
resamp_alg = self.rspAlg_DS,
cliptoextent = cliptoextent,
#clipextent = self.im2shift.box.boxMapYX,
progress = self.progress,
v = self.v,
q = self.q)
path_out = self.path_out,
fmt_out = self.fmt_out,
out_crea_options = self.out_creaOpt,
align_grids = self.align_grids,
match_gsd = self.match_gsd,
out_gsd = self.out_gsd,
target_xyGrid = self.target_xyGrid,
min_points_local_corr = min_points_local_corr,
resamp_alg = self.rspAlg_DS,
cliptoextent = cliptoextent,
#clipextent = self.im2shift.box.boxMapYX,
progress = self.progress,
v = self.v,
q = self.q)
self.deshift_results = DS.correct_shifts()
return self.deshift_results
......
......@@ -52,6 +52,9 @@ class DESHIFTER(object):
default = False
- target_xyGrid(list): a list with an x-grid and a y-grid like [[15,45], [15,45]].
This overrides 'out_gsd', 'align_grids' and 'match_gsd'.
- min_points_local_corr number of valid tie points, below which a global shift correction is performed
instead of a local correction (global X/Y shift is then computed as the mean shift
of the remaining points)(default: 5 tie points)
- resamp_alg(str) the resampling algorithm to be used if neccessary
(valid algorithms: nearest, bilinear, cubic, cubic_spline, lanczos, average, mode,
max, min, med, q1, q3)
......@@ -88,6 +91,7 @@ class DESHIFTER(object):
self.band2process = self.band2process-1 if self.band2process is not None else None # internally handled as band index
self.nodata = kwargs.get('nodata' , self.im2shift.nodata)
self.align_grids = kwargs.get('align_grids' , False)
self.min_points_local_corr = kwargs.get('min_points_local_corr', 5)
self.rspAlg = kwargs.get('resamp_alg' , 'cubic') # TODO accept also integers
self.cliptoextent = kwargs.get('cliptoextent' , False)
self.clipextent = kwargs.get('clipextent' , None)
......@@ -101,18 +105,31 @@ class DESHIFTER(object):
self.shift_prj = self.im2shift.projection
self.shift_gt = list(self.im2shift.geotransform)
# in case of local shift correction and local coreg results contain less points than min_points_local_corr:
# force global correction based on mean X/Y shifts
if 'GCPList' in coreg_results and len(coreg_results['GCPList']) < self.min_points_local_corr:
warnings.warn('Only %s valid tie point(s) could be identified. A local shift correction is therefore not '
'reasonable and could cause artifacts in the output image. The target image is '
'corrected globally with the mean X/Y shift of %.3f/%.3f pixels.' %(len(self.GCPList),
coreg_results['mean_shifts_px']['x'], coreg_results['mean_shifts_px']['y']))
self.GCPList = None
coreg_results['updated map info'] = coreg_results['updated map info means']
# in case of global shift correction -> the updated map info from coreg_results already has the final map info
# BUT: this will updated in correct_shifts() if clipextent is given or warping is needed
if not self.GCPList:
# in case of global de-shifting -> the updated map info from coreg_results already has the final map info
# BUT: this will updated in correct_shifts() if clipextent is given or warping is needed
mapI = coreg_results['updated map info']
self.updated_map_info = mapI if mapI else geotransform2mapinfo(self.shift_gt, self.shift_prj)
self.updated_gt = mapinfo2geotransform(self.updated_map_info) if mapI else self.shift_gt
self.original_map_info = coreg_results['original map info'] # only exists in global de-shifting # FIXME WHY?
self.original_map_info = coreg_results['original map info']
self.updated_projection = self.ref_prj
self.out_grid = self._get_out_grid() # needs self.ref_grid, self.im2shift
self.out_gsd = [abs(self.out_grid[0][1]-self.out_grid[0][0]), abs(self.out_grid[1][1]-self.out_grid[1][0])] # xgsd, ygsd
# assertions
assert self.rspAlg in _dict_rspAlg_rsp_Int.keys(), \
"'%s' is not a supported resampling algorithm." %self.rspAlg
......@@ -122,6 +139,7 @@ class DESHIFTER(object):
'bands' if self.im2shift.bands > 1 else 'band', 'between 1 and ' if self.im2shift.bands > 1 else '',
self.im2shift.bands, self.band2process+1)
# set defaults for general class attributes
self.is_shifted = False # this is not included in COREG.coreg_info
self.is_resampled = False # this is not included in COREG.coreg_info
......
......@@ -99,17 +99,29 @@ class Tie_Point_Grid(object):
self.ref = self.COREG_obj.ref
self.shift = self.COREG_obj.shift
self.XY_points, self.XY_mapPoints = self._get_imXY__mapXY_points(self.grid_res)
self._CoRegPoints_table = None # set by self.CoRegPoints_table
self._GCPList = None # set by self.to_GCPList()
self.kriged = None # set by Raster_using_Kriging()
mean_x_shift_px = property(lambda self:
self.CoRegPoints_table['X_SHIFT_PX'][self.CoRegPoints_table['X_SHIFT_PX']!=self.outFillVal].mean())
mean_y_shift_px = property(lambda self:
self.CoRegPoints_table['Y_SHIFT_PX'][self.CoRegPoints_table['Y_SHIFT_PX']!=self.outFillVal].mean())
mean_x_shift_map = property(lambda self:
self.CoRegPoints_table['X_SHIFT_M'][self.CoRegPoints_table['X_SHIFT_M']!=self.outFillVal].mean())
mean_y_shift_map = property(lambda self:
self.CoRegPoints_table['Y_SHIFT_M'][self.CoRegPoints_table['Y_SHIFT_M']!=self.outFillVal].mean())
@property
def CoRegPoints_table(self):
"""Returns a GeoDataFrame with the columns 'geometry','POINT_ID','X_IM','Y_IM','X_UTM','Y_UTM','X_WIN_SIZE',
'Y_WIN_SIZE','X_SHIFT_PX','Y_SHIFT_PX', 'X_SHIFT_M', 'Y_SHIFT_M', 'ABS_SHIFT' and 'ANGLE' containing all
information containing all the results frm coregistration for all points in the tie points grid.
information containing all the results from coregistration for all points in the tie points grid.
"""
if self._CoRegPoints_table is not None:
return self._CoRegPoints_table
......
......@@ -4,8 +4,8 @@
__author__ = """Daniel Scheffler"""
__email__ = 'daniel.scheffler@gfz-potsdam.de'
__version__ = '0.3.1'
__versionalias__ = '2017-07-04_02'
__version__ = '0.3.4'
__versionalias__ = '2017-07-07_01'
from .CoReg import COREG
......
......@@ -24,7 +24,7 @@ test_requirements = ['coverage']
setup(
name='arosics',
version='0.3.1',
version='0.3.4',
description="An Automated and Robust Open-Source Image Co-Registration Software for Multi-Sensor Satellite Data",
long_description=readme + '\n\n' + history,
author="Daniel Scheffler",
......
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