Commit 55989a50 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

added tie point filtering by measuring SSIM improvement; fastened grid initialization

components.CoReg.imParams:
- added attribute 'xygrid_specs'

components.CoReg.COREG:
- added attributes 'ssim_orig', 'ssim_deshifted', '_ssim_improved' and property 'ssim_improved'
- _get_opt_winpos_winsize(): outlier assertion is now a tracked exception
- _get_image_windows_to_match(): bugfix for passing rspAlg from wrong attribute
- added function _get_deshifted_otherWin()
- added function _validate_ssim_improvement()
- calculate_spatial_shifts() added call for _validate_ssim_improvement()
- added function _get_inverted_coreg_info()

components.CoReg_local.COREG_LOCAL:
- implemented keyword 'tieP_filter_level'
- view_CoRegPoints(): implemented keyword 'hide_filtered'

components.DeShifter.DESHIFTER:
- bugfix for not handling attribute 'band2process' correctly
- added attributes init_args and init_kwargs
- revised _get_out_extent(): bugfix for not handling clip extent correctly; bugfix for wrong call of find_nearest()
- correct_shifts(): bugfix for wrong condition in connection with target_xyGrid keyword

components.Geom_Quality_Grid.Geom_Quality_Grid:
- implemented keyword 'tieP_filter_level'
- get_CoRegPoints_table():
     - exclusion of points outside of overlap polygon is now much faster
     - added new columns to CoRegPoints_table: 'SSIM_BEFORE', 'SSIM_AFTER', 'SSIM_IMPROVED', 'LAST_ERR'

- updated __version__
parent 6504101e
......@@ -9,7 +9,7 @@ from .components import utilities
from .components import geometry
__author__ = 'Daniel Scheffler'
__version__= '2016-11-09_01'
__version__= '2016-11-12_01'
__all__=['COREG',
'COREG_LOCAL',
......
......@@ -37,6 +37,7 @@ from py_tools_ds.ptds.geo.coord_trafo import pixelToMapYX, reproject_shapel
from py_tools_ds.ptds.geo.raster.reproject import warp_ndarray
from py_tools_ds.ptds.geo.map_info import geotransform2mapinfo
from py_tools_ds.ptds.numeric.vector import find_nearest
from py_tools_ds.ptds.similarity.raster import calc_ssim
......@@ -62,13 +63,14 @@ class imParamObj(object):
self.title = os.path.basename(self.GeoArray.filePath) if self.GeoArray.filePath else self.imName
# set params
self.prj = self.GeoArray.projection
self.gt = self.GeoArray.geotransform
self.xgsd = self.GeoArray.xgsd
self.ygsd = self.GeoArray.ygsd
self.rows = self.GeoArray.rows
self.cols = self.GeoArray.cols
self.bands = self.GeoArray.bands
self.prj = self.GeoArray.projection
self.gt = self.GeoArray.geotransform
self.xgsd = self.GeoArray.xgsd
self.ygsd = self.GeoArray.ygsd
self.xygrid_specs = self.GeoArray.xygrid_specs
self.rows = self.GeoArray.rows
self.cols = self.GeoArray.cols
self.bands = self.GeoArray.bands
# validate params
assert self.prj, 'The %s has no projection.' % self.imName
......@@ -77,9 +79,9 @@ class imParamObj(object):
# set band4match
self.band4match = (CoReg_params['r_b4match'] if imID == 'ref' else CoReg_params['s_b4match'])-1
assert self.bands >= self.band4match+1 >= 1, "The %s has %s %s. So its band number to match must be %s%s." \
% (self.imName, self.bands, 'bands' if self.bands > 1 else 'band', 'between 1 and '
if self.bands > 1 else '', self.bands)
assert self.bands >= self.band4match+1 >= 1, "The %s has %s %s. So its band number to match must be %s%s. " \
"Got %s." % (self.imName, self.bands, 'bands' if self.bands > 1 else 'band', 'between 1 and '
if self.bands > 1 else '', self.bands, self.band4match)
# set nodata
if CoReg_params['nodata'][0 if imID == 'ref' else 1] is not None:
......@@ -242,6 +244,9 @@ class COREG(object):
self.vec_length_map = None
self.vec_angle_deg = None
self.updated_map_info = None # set by self.get_updated_map_info()
self.ssim_orig = None # set by self._validate_ssim_improvement()
self.ssim_deshifted = None # set by self._validate_ssim_improvement()
self._ssim_improved = None # private attribute to be filled by self.ssim_improved
self.tracked_errors = [] # expanded each time an error occurs
self.success = None # default
......@@ -394,8 +399,14 @@ class COREG(object):
wp = (wp[0] if wp[0] else overlap_center_pos_x[0]), (wp[1] if wp[1] else overlap_center_pos_y[0])
# validate window position
assert self.overlap_poly.contains(Point(wp)), 'The provided window position %s/%s is outside of the overlap ' \
'area of the two input images. Check the coordinates.' %wp
if not self.overlap_poly.contains(Point(wp)):
self.success=False
self.tracked_errors.append(ValueError('The provided window position %s/%s is outside of the overlap ' \
'area of the two input images. Check the coordinates.' %wp))
if not self.ignErr:
raise self.tracked_errors[-1]
# assert self.overlap_poly.contains(Point(wp)), 'The provided window position %s/%s is outside of the overlap ' \
# 'area of the two input images. Check the coordinates.' %wp
#for im in [self.ref, self.shift]:
# imX, imY = mapXY2imXY(wp, im.gt)
#if self.ignErr:
......@@ -522,7 +533,7 @@ class COREG(object):
self.matchWin.imParams.prj,
out_gsd = (self.imfft_gsd, self.imfft_gsd),
out_bounds = ([tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax]),
rspAlg = self.rspAlg_calc,
rspAlg = _dict_rspAlg_rsp_Int[self.rspAlg_calc],
in_nodata = self.otherWin.imParams.nodata,
CPUs = None if self.mp else 1,
progress = False) [0]
......@@ -751,78 +762,110 @@ class COREG(object):
return x_intshift+x_subshift, y_intshift+y_subshift
def _validate_ssim_improvement(self):
def _get_deshifted_otherWin(self):
"""Returns a de-shifted version of self.otherWin as a GeoArray instance.The output dimensions and geographic
bounds are equal to those of self.matchWin and geometric shifts are corrected according to the previously
computed X/Y shifts within the matching window. This allows direct application of algorithms e.g. measuring
image similarity.
The image subset that is resampled in this function is always the same that has been resampled during
computation of geometric shifts (usually the image with the higher geometric resolution).
:returns: GeoArray instance of de-shifted self.otherWin
"""
# shift vectors have been calculated to fit target image onto reference image
# -> so the shift vectors have to be inverted if shifts are applied to reference image
coreg_info = self._get_inverted_coreg_info() if self.otherWin.imParams.imName=='reference image' else \
self.coreg_info
ds_results = DESHIFTER(self.otherWin.imParams.GeoArray, coreg_info,
band2process = self.otherWin.imParams.band4match+1,
clipextent = list(np.array(self.matchWin.boundsMap)[[0,2,1,3]]),
target_xyGrid = self.matchWin.imParams.xygrid_specs,
q = True
).correct_shifts()
return ds_results['GeoArray_shifted']
def _validate_ssim_improvement(self, v=False):
"""Computes mean structural similarity index between reference and target image before and after correction
of geometric shifts..
:param v: <bool> verbose mode: shows images of the matchWin, otherWin and shifted version of otherWin
:return: <tuple> SSIM before an after shift correction
"""
assert self.success is not None,\
'Calculate geometric shifts first before trying to measure image similarity improvement!'
assert self.success in [True, None],\
'Since calculation of geometric shifts failed, no image similarity improvement can be measured.'
# get image dynamic range
dr = max(self.ref.win.data.max(), self.shift.win.data.max()) - \
min(self.ref.win.data.max(), self.shift.win.data.max())
# compute ssim BEFORE shift correction
from py_tools_ds.ptds.similarity.raster import calc_ssim
ssim_before = calc_ssim(self.matchWin.data, self.otherWin.data, dynamic_range=dr)
print('SSIM before', ssim_before)
#ws = int(self.matchWin.imDimsYX[1]), int(self.matchWin.imDimsYX[0])
# get shifted GeoArray in the reference image grid
shifted_geoArr = copy(self.shift.GeoArray)
geotransform = list(shifted_geoArr.gt)
geotransform[0] += self.x_shift_map
geotransform[3] += self.y_shift_map
shifted_geoArr.gt = geotransform
tgt_xmin, tgt_xmax, tgt_ymin, tgt_ymax = self.ref.win.boundsMap
arr2warp = shifted_geoArr[:,:,self.shift.band4match] if shifted_geoArr.ndim==3 else shifted_geoArr[:] # FIXME dont read complete band
from py_tools_ds.ptds.io.raster.GeoArray import get_array_at_mapPos
sub_arr, sub_gt, sub_prj = get_array_at_mapPos(arr2warp, shifted_geoArr.gt, shifted_geoArr.prj, self.ref.prj,
mapBounds=(tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax),
mapBounds_prj=self.ref.prj,
out_gsd=(5,5), # FIXME stimmt das?
rspAlg='cubic')
# # otherWin per subset-read einlesen
# rS, rE, cS, cE = GEO.get_GeoArrayPosition_from_boxImYX(self.otherWin.boxImYX)
# assert np.array_equal(np.abs(np.array([rS, rE, cS, cE])), np.array([rS, rE, cS, cE])), \
# 'Got negative values in gdalReadInputs for %s.' % self.otherWin.imParams.imName
# data2warp = shifted_geoArr[rS:rE, cS:cE, self.otherWin.imParams.band4match]
#
#
# if self.v:
# print('Original matching windows:')
# ref_data, shift_data = (self.matchWin.data, self.otherWin.data) if self.grid2use=='ref' else \
# (self.otherWin.data, self.matchWin.data)
# PLT.subplot_imshow([ref_data, shift_data],[self.ref.title,self.shift.title], grid=True)
#
# otherWin_subgt = GEO.get_subset_GeoTransform(self.otherWin.gt, self.otherWin.boxImYX)
#
# # resample otherWin.data to the resolution of matchWin AND make sure the pixel edges are identical
# # (in order to make each image show the same window with the same coordinates)
# # TODO replace cubic resampling by PSF resampling - average resampling leads to sinus like distortions in the fft image that make a precise coregistration impossible. Thats why there is currently no way around cubic resampling.
# #tgt_xmin,tgt_xmax,tgt_ymin,tgt_ymax = self.matchWin.boundsMap
# sub_arr = warp_ndarray(data2warp,
# otherWin_subgt,
# self.shift.prj,
# self.ref.prj,
# out_gsd = (self.imfft_gsd, self.imfft_gsd),
# out_bounds = ([tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax]),
# rspAlg = self.rspAlg_calc,
# in_nodata = self.shift.nodata,
# CPUs = None if self.mp else 1,
# progress = False) [0]
#out_arr, out_gt, out_prj = \
# warp_ndarray(arr, arr_gt, arr_prj, out_prj=out_prj, out_bounds=mapBounds, out_bounds_prj=mapBounds_prj,
# in_nodata=fillVal, out_nodata=fillVal, rspAlg=rspAlg, out_gsd=out_gsd)
print(sub_arr.shape)
#GeoArray(sub_arr).show(figsize=(15,15))
ssim_after = calc_ssim(sub_arr, self.otherWin.data)
print('SSIM after', ssim_after)
dr = max(self.matchWin.data.max(), self.otherWin.data.max()) - \
min(self.matchWin.data.min(), self.otherWin.data.min())
# compute SSIM BEFORE shift correction
self.ssim_orig = calc_ssim(self.matchWin.data, self.otherWin.data, dynamic_range=dr, gaussian_weights=True)
# compute SSIM AFTER shift correction
## resample otherWin while correcting detected shifts and match geographic bounds of matchWin
otherWin_deshift_geoArr = self._get_deshifted_otherWin()
## get the corresponding matchWin data
matchWinData = self.matchWin.data
## check if shapes of two images are equal (due to bug (?), in some cases otherWin_deshift_geoArr does not have
## the exact same dimensions as self.matchWin -> maybe bounds are handled differently by gdal.Warp)
if not self.matchWin.data.shape == otherWin_deshift_geoArr.shape:
matchWinData, matchWinGt, matchWinPrj = self.matchWin.imParams.GeoArray.get_mapPos(
list(np.array(self.matchWin.boundsMap)[[0, 2, 1, 3]]), self.matchWin.imParams.prj, rspAlg='cubic',
band2get=self.matchWin.imParams.band4match)
self.ssim_deshifted = calc_ssim(otherWin_deshift_geoArr[:], matchWinData, dynamic_range=dr, gaussian_weights=True)
if v:
GeoArray(matchWinData).show()
GeoArray(self.otherWin.data).show()
otherWin_deshift_geoArr.show()
if not self.q:
print('Image similarity within the matching window (SSIM before/after correction): %.4f => %.4f'
% (self.ssim_orig, self.ssim_deshifted))
self.ssim_improved = self.ssim_orig < self.ssim_deshifted
# write win data to disk
#outDir = '/home/gfz-fe/scheffler/temp/ssim_debugging/'
#GeoArray(matchWinData, matchWinGt, matchWinPrj).save(outDir+'matchWinData.bsq')
#otherWinGt = (self.otherWin.boundsMap[0], self.matchWin.imParams.xgsd, 0, self.otherWin.boundsMap[3], 0, -self.matchWin.imParams.ygsd)
#GeoArray(self.otherWin.data, therWinGt, self.otherWin.imParams.prj).save(outDir+'otherWin.data.bsq')
# otherWin_deshift_geoArr.save(outDir+''shifted.bsq')
return self.ssim_orig, self.ssim_deshifted
@property
def ssim_improved(self):
"""Returns True if image similarity within the matching window has been improved by correcting the previously
computed geometric shifts."""
if self.success is True:
if self._ssim_improved is None:
ssim_orig, ssim_deshifted = self._validate_ssim_improvement()
self._ssim_improved = ssim_orig <= ssim_deshifted
return self._ssim_improved
@ssim_improved.setter
def ssim_improved(self, has_improved):
self._ssim_improved = has_improved
def calculate_spatial_shifts(self):
......@@ -933,9 +976,13 @@ class COREG(object):
print('Calculated absolute shift vector length in map units: %s' %self.vec_length_map)
print('Calculated angle of shift vector in degrees from North: %s' %self.vec_angle_deg)
if self.x_shift_px or self.y_shift_px:
self._get_updated_map_info()
# set self.ssim_before and ssim_after
self._validate_ssim_improvement()
warnings.simplefilter('default')
return 'success'
......@@ -971,6 +1018,28 @@ class COREG(object):
return self.coreg_info
def _get_inverted_coreg_info(self):
"""Returns an inverted dictionary of coreg_info that can be passed to DESHIFTER in order to fit the REFERENCE
image onto the TARGET image."""
inv_coreg_info = copy(self.coreg_info)
inv_coreg_info['corrected_shifts_px']['x'] *= -1
inv_coreg_info['corrected_shifts_px']['y'] *= -1
inv_coreg_info['corrected_shifts_map']['x'] *= -1
inv_coreg_info['corrected_shifts_map']['y'] *= -1
inv_coreg_info['original map info'] = geotransform2mapinfo(self.ref.gt, self.ref.prj)
inv_coreg_info['reference geotransform'] = self.shift.gt
inv_coreg_info['reference grid'] = self.shift.xygrid_specs
inv_coreg_info['reference projection'] = self.shift.prj
if inv_coreg_info['updated map info']:
updated_map_info = copy(inv_coreg_info['original map info'] )
updated_map_info[3] = str(float(inv_coreg_info['original map info'][3]) - self.x_shift_map)
updated_map_info[4] = str(float(inv_coreg_info['original map info'][4]) - self.y_shift_map)
inv_coreg_info['updated map info'] = updated_map_info
return inv_coreg_info
def correct_shifts(self):
DS = DESHIFTER(self.shift.GeoArray, self.coreg_info,
path_out = self.path_out,
......
......@@ -26,10 +26,10 @@ class COREG_LOCAL(object):
def __init__(self, im_ref, im_tgt, grid_res, window_size=(256,256), path_out=None, fmt_out='ENVI',
out_crea_options=None, projectDir=None, r_b4match=1, s_b4match=1, max_iter=5, max_shift=5,
align_grids=True, match_gsd=False, out_gsd=None, target_xyGrid=None, resamp_alg_deshift='cubic',
resamp_alg_calc='cubic', footprint_poly_ref=None, footprint_poly_tgt=None, data_corners_ref=None,
data_corners_tgt=None, outFillVal=-9999, nodata=(None, None), calc_corners=True, binary_ws=True,
CPUs=None, progress=True, v=False, q=False, ignore_errors=False):
tieP_filter_level=1, align_grids=True, match_gsd=False, out_gsd=None, target_xyGrid=None,
resamp_alg_deshift='cubic', resamp_alg_calc='cubic', footprint_poly_ref=None, footprint_poly_tgt=None,
data_corners_ref=None, data_corners_tgt=None, outFillVal=-9999, nodata=(None, None), calc_corners=True,
binary_ws=True, CPUs=None, progress=True, v=False, q=False, ignore_errors=False):
"""Applies the algorithm to detect spatial shifts to the whole overlap area of the input images. Spatial shifts
are calculated for each point in grid of which the parameters can be adjusted using keyword arguments. Shift
......@@ -53,6 +53,11 @@ class COREG_LOCAL(object):
:param s_b4match(int): band of shift image to be used for matching (starts with 1; default: 1)
:param max_iter(int): maximum number of iterations for matching (default: 5)
:param max_shift(int): maximum shift distance in reference image pixel units (default: 5 px)
:param tieP_filter_level(int): filter tie points used for shift correction in different levels:
- Level 0: no tie point filtering
- Level 1: SSIM filtering - filters all tie points out where shift
correction does not increase image similarity within matching window
(measured by mean structural similarity index)
:param out_gsd (float): output pixel size in units of the reference coordinate system (default = pixel
size of the input array), given values are overridden by match_gsd=True
:param align_grids (bool): True: align the input coordinate grid to the reference (does not affect the
......@@ -120,6 +125,7 @@ class COREG_LOCAL(object):
self.window_size = window_size
self.max_shift = max_shift
self.max_iter = max_iter
self.tieP_filter_level = tieP_filter_level
self.align_grids = align_grids
self.match_gsd = match_gsd
self.out_gsd = out_gsd
......@@ -191,9 +197,15 @@ class COREG_LOCAL(object):
if self._quality_grid:
return self._quality_grid
else:
self._quality_grid = Geom_Quality_Grid(self.COREG_obj, self.grid_res, outFillVal=self.outFillVal,
resamp_alg_calc=self.rspAlg_calc, dir_out=self.projectDir,
CPUs=self.CPUs, progress=self.progress, v=self.v, q=self.q)
self._quality_grid = Geom_Quality_Grid(self.COREG_obj, self.grid_res,
outFillVal = self.outFillVal,
resamp_alg_calc = self.rspAlg_calc,
tieP_filter_level = self.tieP_filter_level,
dir_out = self.projectDir,
CPUs = self.CPUs,
progress = self.progress,
v = self.v,
q = self.q)
if self.v:
self.view_CoRegPoints(figsize=(10,10))
return self._quality_grid
......@@ -223,7 +235,7 @@ class COREG_LOCAL(object):
def view_CoRegPoints(self, attribute2plot='ABS_SHIFT', cmap=None, exclude_fillVals=True, backgroundIm='tgt',
figsize=None, savefigPath='', savefigDPI=96):
hide_filtered=True, figsize=None, savefigPath='', savefigDPI=96):
"""Shows a map of the calculated quality grid with the target image as background.
:param attribute2plot: <str> the attribute of the quality grid to be shown (default: 'ABS_SHIFT')
......@@ -232,6 +244,7 @@ class COREG_LOCAL(object):
:param exclude_fillVals: <bool> whether to exclude those points of the grid where spatial shift detection failed
:param backgroundIm: <str> whether to use the target or the reference image as map background. Possible
options are 'ref' and 'tgt' (default: 'tgt')
:param hide_filtered: <bool> hide all points that have been filtered out according to tie point filter level
:param figsize: <tuple> size of the figure to be viewed, e.g. (10,10)
:param savefigPath:
:param savefigDPI:
......@@ -248,8 +261,8 @@ class COREG_LOCAL(object):
# transform all points of quality grid to LonLat
GDF = self.CoRegPoints_table.loc\
[self.CoRegPoints_table.X_SHIFT_M != self.outFillVal, ['geometry', attribute2plot]].copy() \
if exclude_fillVals else self.CoRegPoints_table.loc[:, ['geometry', attribute2plot]]
[self.CoRegPoints_table.X_SHIFT_M != self.outFillVal, ['geometry', attribute2plot, 'SSIM_IMPROVED']].copy() \
if exclude_fillVals else self.CoRegPoints_table.loc[:, ['geometry', attribute2plot, 'SSIM_IMPROVED']]
# get LonLat coordinates for all points
get_LonLat = lambda X, Y: transform_any_prj(self.im2shift.projection, 4326, X, Y)
......@@ -268,6 +281,16 @@ class COREG_LOCAL(object):
GDF['plt_XY'] = list(GDF['LonLat'].map(lambda ll: map2show(*ll)))
GDF['plt_X'] = list(GDF['plt_XY'].map(lambda XY: XY[0]))
GDF['plt_Y'] = list(GDF['plt_XY'].map(lambda XY: XY[1]))
if not hide_filtered and self.tieP_filter_level>0:
# flag SSIM filtered points
GDF_filt = GDF[GDF.SSIM_IMPROVED==False].copy()
plt.scatter(GDF_filt['plt_X'], GDF_filt['plt_Y'], c='r', marker='o' if len(GDF) < 10000 else '.', s=150, alpha=1.0)
if hide_filtered:
GDF = GDF[GDF.SSIM_IMPROVED==True].copy()
# plot all points on top
vmin, vmax = np.percentile(GDF[attribute2plot], 0), np.percentile(GDF[attribute2plot], 95)
points = plt.scatter(GDF['plt_X'],GDF['plt_Y'], c=GDF[attribute2plot],
cmap=palette, marker='o' if len(GDF)<10000 else '.', s=50, alpha=1.0,
......
......@@ -63,6 +63,10 @@ class DESHIFTER(object):
- q(bool): quiet mode (default: False)
"""
# TODO validate kwargs for typos
self.init_args = dict([x for x in locals().items() if x[0] != "self" and not x[0].startswith('__')])
self.init_kwargs = self.init_args['kwargs']
# unpack args
self.im2shift = im2shift if isinstance(im2shift, GeoArray) else GeoArray(im2shift)
self.GCPList = coreg_results['GCPList'] if 'GCPList' in coreg_results else None
......@@ -75,6 +79,7 @@ class DESHIFTER(object):
self.fmt_out = kwargs.get('fmt_out' , 'ENVI')
self.out_creaOpt = kwargs.get('out_crea_options', [])
self.band2process = kwargs.get('band2process' , None) # starts with 1 # FIXME warum?
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.rspAlg = kwargs.get('resamp_alg' , 'cubic')
......@@ -89,6 +94,7 @@ class DESHIFTER(object):
self.im2shift.q = self.q
self.shift_prj = self.im2shift.projection
self.shift_gt = list(self.im2shift.geotransform)
if not self.GCPList:
mapI = coreg_results['updated map info']
self.updated_map_info = mapI if mapI else geotransform2mapinfo(self.shift_gt, self.shift_prj)
......@@ -96,12 +102,17 @@ class DESHIFTER(object):
self.updated_gt = mapinfo2geotransform(self.updated_map_info) if mapI else self.shift_gt
self.updated_projection = self.ref_prj
self.out_grid = self._get_out_grid(kwargs) # needs self.ref_grid, self.im2shift
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
if self.band2process is not None:
assert self.im2shift.bands-1 >= self.band2process >= 0, "The %s '%s' has %s %s. So 'band2process' must be " \
"%s%s. Got %s." % (self.im2shift.__class__.__name__, self.im2shift.basename, self.im2shift.bands,
'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
......@@ -111,11 +122,11 @@ class DESHIFTER(object):
self.GeoArray_shifted = None # set by self.correct_shifts
def _get_out_grid(self, init_kwargs):
def _get_out_grid(self):
# parse given params
out_gsd = init_kwargs.get('out_gsd' , None)
match_gsd = init_kwargs.get('match_gsd' , False)
out_grid = init_kwargs.get('target_xyGrid', None)
out_gsd = self.init_kwargs.get('out_gsd' , None)
match_gsd = self.init_kwargs.get('match_gsd' , False)
out_grid = self.init_kwargs.get('target_xyGrid', None)
# assertions
assert out_grid is None or (isinstance(out_grid,(list, tuple)) and len(out_grid)==2)
......@@ -159,35 +170,37 @@ class DESHIFTER(object):
return get_grid(self.im2shift.geotransform, self.im2shift.xgsd, self.im2shift.ygsd)
@staticmethod
def _grids_alignable(in_xgsd, in_ygsd, out_xgsd, out_ygsd):
def _grids_alignable(self,in_xgsd, in_ygsd, out_xgsd, out_ygsd):
is_alignable = lambda gsd1, gsd2: max(gsd1, gsd2) % min(gsd1, gsd2) == 0 # checks if pixel sizes are divisible
if not is_alignable(in_xgsd, out_xgsd) or not is_alignable(in_ygsd, out_ygsd):
warnings.warn("\nThe targeted output coordinate grid is not alignable with the image to be shifted because "
"their pixel sizes are not exact multiples of each other (input [X/Y]: "
"%s %s; output [X/Y]: %s %s). Therefore the targeted output grid is "
"chosen for the resampled output image. If you don´t like that you can use the '-out_gsd' "
"parameter to set an appropriate output pixel size.\n"
% (in_xgsd, in_ygsd, out_xgsd, out_ygsd))
if not self.q:
warnings.warn("\nThe coordinate grid of the image to be shifted cannot be aligned to the reference "
"grid because their pixel sizes are not exact multiples of each other (input [X/Y]: "
"%s/%s; output [X/Y]: %s/%s). Therefore the original grid is chosen for the resampled "
"output image. If you don´t like that you can use the '-out_gsd' parameter to set an "
"appropriate output pixel size.\n" % (in_xgsd, in_ygsd, out_xgsd, out_ygsd))
return False
else:
return True
def _get_out_extent(self):
if self.cliptoextent and self.clipextent is None:
self.clipextent = self.im2shift.footprint_poly.envelope.bounds
else:
xmin, xmax, ymin, ymax = self.im2shift.box.boundsMap
self.clipextent = xmin, ymin, xmax, ymax
if self.clipextent is None:
# no clip extent has been given
if self.cliptoextent:
# use actual image corners as clip extent
self.clipextent = self.im2shift.footprint_poly.envelope.bounds
else:
# use outer bounds of the image as clip extent
xmin, xmax, ymin, ymax = self.im2shift.box.boundsMap
self.clipextent = xmin, ymin, xmax, ymax
# snap clipextent to output grid (in case of odd input coords the output coords are moved INSIDE the input array)
xmin, ymin, xmax, ymax = self.clipextent
xmin = find_nearest(self.out_grid[0], xmin, roundAlg='on' , extrapolate=True)
ymin = find_nearest(self.out_grid[1], ymin, roundAlg='on' , extrapolate=True)
xmax = find_nearest(self.out_grid[0], xmax, roundAlg='off', extrapolate=True)
ymax = find_nearest(self.out_grid[0], ymax, roundAlg='off', extrapolate=True)
ymax = find_nearest(self.out_grid[1], ymax, roundAlg='off', extrapolate=True)
return xmin, ymin, xmax, ymax
......@@ -200,7 +213,8 @@ class DESHIFTER(object):
t_start = time.time()
equal_prj = prj_equal(self.ref_prj,self.shift_prj)
if equal_prj and is_coord_grid_equal(self.shift_gt, *self.out_grid) and not self.align_grids and not self.GCPList:
if equal_prj and is_coord_grid_equal(self.shift_gt, *self.out_grid) and not self.align_grids and \
not self.GCPList and not self.init_kwargs.get('target_xyGrid',None):
# FIXME buggy condition:
# reconstructable with correct_spatial_shifts from GMS
#DS = DESHIFTER(geoArr, self.coreg_info,
......@@ -217,11 +231,12 @@ class DESHIFTER(object):
# clip with target extent
self.arr_shifted, self.updated_gt, self.updated_projection = \
shifted_geoArr.get_mapPos((xmin,ymin,xmax,ymax), self.shift_prj, fillVal=self.nodata)
shifted_geoArr.get_mapPos((xmin,ymin,xmax,ymax), self.shift_prj, fillVal=self.nodata,
band2get=self.band2process)
self.updated_map_info = geotransform2mapinfo(self.updated_gt, self.updated_projection)
else:
# array keeps the same; updated gt and prj are taken from coreg_info
self.arr_shifted = self.im2shift[:]
self.arr_shifted = self.im2shift[:,:,self.band2process]
self.GeoArray_shifted = GeoArray(self.arr_shifted, tuple(self.shift_gt), self.updated_projection)
if self.path_out:
......@@ -270,7 +285,9 @@ class DESHIFTER(object):
#
# # TO DO implement output writer
in_arr = self.im2shift[self.band2process] if self.band2process else self.im2shift[:]
# FIXME avoid reading the whole band if clip_extent is passed
in_arr = self.im2shift[:,:,self.band2process] if self.band2process is not None and self.im2shift.ndim==3\
else self.im2shift[:]
if not self.GCPList:
# apply XY-shifts to shift_gt
......@@ -290,7 +307,7 @@ class DESHIFTER(object):
#warpOptions= ['-refine_gcps 500 1.9'],
#options = '-wm 10000',# -order 3',
#options = ['-order 3'],
# options = ['GDAL_CACHEMAX 800 '],
#options = ['GDAL_CACHEMAX 800 '],
#warpMemoryLimit=125829120, # 120MB
CPUs = self.CPUs,
progress = self.progress,
......
......@@ -28,13 +28,14 @@ from py_tools_ds.ptds.processing.progress_mon import ProgressBar
global_shared_imref = None
global_shared_im2shift = None
overlap_poly =None
class Geom_Quality_Grid(object):
"""See help(Geom_Quality_Grid) for documentation!"""
def __init__(self, COREG_obj, grid_res, outFillVal=-9999, resamp_alg_calc='cubic', dir_out=None, CPUs=None,
progress=True, v=False, q=False):
def __init__(self, COREG_obj, grid_res, outFillVal=-9999, resamp_alg_calc='cubic', tieP_filter_level=1,
dir_out=None, CPUs=None, progress=True, v=False, q=False):
"""Applies the algorithm to detect spatial shifts to the whole overlap area of the input images. Spatial shifts
are calculated for each point in grid of which the parameters can be adjusted using keyword arguments. Shift
......@@ -50,6 +51,11 @@ class Geom_Quality_Grid(object):
(valid algorithms: nearest, bilinear, cubic, cubic_spline, lanczos, average, mode,
max, min, med, q1, q3)
default: cubic (highly recommended)
:param tieP_filter_level(int): filter tie points used for shift correction in different levels:
- Level 0: no tie point filtering
- Level 1: SSIM filtering - filters all tie points out where shift
correction does not increase image similarity within matching window
(measured by mean structural similarity index)
:param dir_out(str): output directory to be used for all outputs if nothing else is given
to the individual methods
:param CPUs(int): number of CPUs to use during calculation of geometric quality grid
......@@ -61,18 +67,19 @@ class Geom_Quality_Grid(object):
if not isinstance(COREG_obj, COREG): raise ValueError("'COREG_obj' must be an instance of COREG class.")
self.COREG_obj = COREG_obj
self.grid_res = grid_res