Commit 596fd660 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'dev'

parents 184ef069 55989a50
......@@ -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',
......
This diff is collapsed.
......@@ -26,9 +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,
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
......@@ -52,6 +53,29 @@ 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
output pixel size as long as input and output pixel sizes are compatible
(5:30 or 10:30 but not 4:30), default = True
:param match_gsd (bool): True: match the input pixel size to the reference pixel size,
default = False
:param target_xyGrid(list): a list with a target x-grid and a target y-grid like [[15,45], [15,45]]
This overrides 'out_gsd', 'align_grids' and 'match_gsd'.
:param resamp_alg_deshift(str) the resampling algorithm to be used for shift correction (if neccessary)
valid algorithms: nearest, bilinear, cubic, cubic_spline, lanczos, average, mode,
max, min, med, q1, q3
default: cubic
:param resamp_alg_calc(str) the resampling algorithm to be used for all warping processes during calculation
of spatial shifts
(valid algorithms: nearest, bilinear, cubic, cubic_spline, lanczos, average, mode,
max, min, med, q1, q3)
default: cubic (highly recommended)
:param footprint_poly_ref(str): footprint polygon of the reference image (WKT string or shapely.geometry.Polygon),
e.g. 'POLYGON ((299999 6000000, 299999 5890200, 409799 5890200, 409799 6000000,
299999 6000000))'
......@@ -76,9 +100,11 @@ class COREG_LOCAL(object):
:param q(bool): quiet mode (default: False)
:param ignore_errors(bool): Useful for batch processing. (default: False)
"""
# TODO outgsd, matchgsd
# assertions
assert fmt_out
assert fmt_out, "'%s' is not a valid GDAL driver code." %fmt_out
if match_gsd and out_gsd: warnings.warn("'-out_gsd' is ignored because '-match_gsd' is set.\n")
if out_gsd: assert isinstance(out_gsd, list) and len(out_gsd) == 2, 'out_gsd must be a list with two values.'
self.params = dict([x for x in locals().items() if x[0] != "self" and not x[0].startswith('__')])
......@@ -99,6 +125,13 @@ 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
self.target_xyGrid = target_xyGrid
self.rspAlg_DS = resamp_alg_deshift
self.rspAlg_calc = resamp_alg_calc
self.calc_corners = calc_corners
self.nodata = nodata
self.outFillVal = outFillVal
......@@ -122,6 +155,7 @@ class COREG_LOCAL(object):
footprint_poly_tgt = footprint_poly_tgt,
data_corners_ref = data_corners_ref,
data_corners_tgt = data_corners_tgt,
resamp_alg_calc = self.rspAlg_calc,
calc_corners = calc_corners,
r_b4match = r_b4match,
s_b4match = s_b4match,
......@@ -163,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,
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
......@@ -195,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')
......@@ -204,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:
......@@ -220,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)
......@@ -240,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,
......@@ -337,8 +388,11 @@ class COREG_LOCAL(object):
path_out = self.path_out,
fmt_out = self.fmt_out,
out_crea_options = self.out_creaOpt,
out_gsd = (self.im2shift.xgsd,self.im2shift.ygsd),
align_grids = True,
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,
......
......@@ -48,7 +48,8 @@ class DESHIFTER(object):
(5:30 or 10:30 but not 4:30), default = False
- match_gsd (bool): True: match the input pixel size to the reference pixel size,
default = False
- target_xyGrid(list): a list with an x-grid and a y-grid like [[15,45], [15,45]]
- 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'.
- 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)
......@@ -62,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
......@@ -74,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')
......@@ -88,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)
......@@ -95,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
......@@ -110,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)
......@@ -158,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
......@@ -199,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,
......@@ -216,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:
......@@ -269,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
......@@ -289,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,12 +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, 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
......@@ -44,6 +46,16 @@ class Geom_Quality_Grid(object):
:param grid_res: grid resolution in pixels of the target image
:param outFillVal(int): if given the generated geometric quality grid is filled with this value in case
no match could be found during co-registration (default: -9999)
:param resamp_alg_calc(str) the resampling algorithm to be used for all warping processes during calculation
of spatial shifts
(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
......@@ -55,17 +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
self.dir_out = dir_out
self.outFillVal = outFillVal
self.CPUs = CPUs
self.v = v
self.q = q if not v else False # overridden by v
self.progress = progress if not q else False # overridden by q
self.COREG_obj = COREG_obj
self.grid_res = grid_res
self.outFillVal = outFillVal
self.rspAlg_calc = resamp_alg_calc
self.tieP_filter_level = tieP_filter_level
self.dir_out = dir_out
self.CPUs = CPUs
self.v = v
self.q = q if not v else False # overridden by v
self.progress = progress if not q else False # overridden by q
self.ref = self.COREG_obj.ref .GeoArray
self.shift = self.COREG_obj.shift.GeoArray
self.ref = self.COREG_obj.ref .GeoArray
self.shift = self.COREG_obj.shift.GeoArray
self.XY_points, self.XY_mapPoints = self._get_imXY__mapXY_points(self.grid_res)
self._CoRegPoints_table = None # set by self.CoRegPoints_table
......@@ -143,13 +157,18 @@ class Geom_Quality_Grid(object):
assert global_shared_im2shift is not None
CR = COREG(global_shared_imref, global_shared_im2shift, multiproc=False, **coreg_kwargs)
CR.calculate_spatial_shifts()
CR_res = [int(CR.matchWin.imDimsYX[1]), int(CR.matchWin.imDimsYX[0]),
CR.x_shift_px, CR.y_shift_px, CR.x_shift_map, CR.y_shift_map,
CR.vec_length_map, CR.vec_angle_deg]
last_err = CR.tracked_errors[-1] if CR.tracked_errors else None
win_sz_y, win_sz_x = CR.matchWin.imDimsYX if CR.matchWin else (None, None)
CR_res = [win_sz_x, win_sz_y,
CR.x_shift_px, CR.y_shift_px, CR.x_shift_map, CR.y_shift_map,
CR.vec_length_map, CR.vec_angle_deg, CR.ssim_orig, CR.ssim_deshifted, CR.ssim_improved, last_err]
return [pointID]+CR_res
@staticmethod
def is_within(point):
return point.within(overlap_poly)
def get_CoRegPoints_table(self, exclude_outliers=1):
assert self.XY_points is not None and self.XY_mapPoints is not None
......@@ -166,7 +185,7 @@ class Geom_Quality_Grid(object):
# self.path_im2shift = tgt_pathTmp
#ref_ds=tgt_ds=None
XYarr2PointGeom = np.vectorize(lambda X,Y: Point(X,Y), otypes=[Point])
XYarr2PointGeom = np.vectorize(lambda X,Y: Point(X,Y), otypes=[Point]) # FIXME replace that with coord image reprojection
geomPoints = np.array(XYarr2PointGeom(self.XY_mapPoints[:,0],self.XY_mapPoints[:,1]))
if isProjectedOrGeographic(self.COREG_obj.shift.prj)=='geographic':
......@@ -174,7 +193,7 @@ class Geom_Quality_Grid(object):
elif isProjectedOrGeographic(self.COREG_obj.shift.prj)=='projected':
UTMzone = abs(get_UTMzone(prj=self.COREG_obj.shift.prj))
south = get_UTMzone(prj=self.COREG_obj.shift.prj)<0
crs = dict(ellps='WGS84', datum='WGS84', proj='utm', zone=UTMzone,south=south,units='m', no_defs=True)
crs = dict(ellps='WGS84', datum='WGS84', proj='utm', zone=UTMzone, south=south, units='m', no_defs=True)
if not south: del crs['south']
else:
crs = None
......@@ -187,7 +206,7 @@ class Geom_Quality_Grid(object):
GDF.loc[:,['X_UTM','Y_UTM']] = self.XY_mapPoints
if exclude_outliers:
GDF = GDF[GDF['geometry'].within(self.COREG_obj.overlap_poly)]
GDF = GDF[GDF['geometry'].within(self.COREG_obj.overlap_poly.simplify(tolerance=15))]
assert not GDF.empty, 'No coregistration point could be placed within the overlap area. Check yout input data!' # FIXME track that
......@@ -212,6 +231,7 @@ class Geom_Quality_Grid(object):
'pointID' : pID,
'wp' : wp,
'ws' : self.COREG_obj.win_size_XY,
'resamp_alg_calc' : self.rspAlg_calc,
'footprint_poly_ref' : self.COREG_obj.ref.poly,
'footprint_poly_tgt' : self.COREG_obj.shift.poly,
'r_b4match' : self.COREG_obj.ref.band4match+1, # band4match is internally saved as index, starting from 0
......@@ -229,9 +249,9 @@ class Geom_Quality_Grid(object):
# run co-registration for whole grid
if self.CPUs is None or self.CPUs>1:
if not self.q:
print("Calculating geometric quality grid (%s points) in mode 'multiprocessing'..." %len(GDF))
cpus = self.CPUs if self.CPUs is not None else multiprocessing.cpu_count()
print("Calculating geometric quality grid (%s points) using %s CPU cores..." %(len(GDF), cpus))
t0 = time.time()
with multiprocessing.Pool(self.CPUs) as pool:
if self.q or not self.progress:
results = pool.map(self._get_spatial_shifts, list_coreg_kwargs)
......@@ -248,8 +268,8 @@ class Geom_Quality_Grid(object):
break
else:
if not self.q:
print("Calculating geometric quality grid (%s points) in mode 'singleprocessing'..." %len(GDF))
results = np.empty((len(geomPoints),9))
print("Calculating geometric quality grid (%s points) 1 CPU core..." %len(GDF))
results = np.empty((len(geomPoints),13), np.object)
bar = ProgressBar(prefix='\tprogress:')
for i,coreg_kwargs in enumerate(list_coreg_kwargs):
if self.progress:
......@@ -258,9 +278,10 @@ class Geom_Quality_Grid(object):
# FIXME in some cases the code hangs here ==> x
# merge results with GDF
records = GeoDataFrame(np.array(results, np.object), columns=['POINT_ID', 'X_WIN_SIZE', 'Y_WIN_SIZE',
'X_SHIFT_PX','Y_SHIFT_PX', 'X_SHIFT_M',
'Y_SHIFT_M', 'ABS_SHIFT', 'ANGLE'])
records = GeoDataFrame(np.array(results, np.object),
columns=['POINT_ID', 'X_WIN_SIZE', 'Y_WIN_SIZE', 'X_SHIFT_PX','Y_SHIFT_PX', 'X_SHIFT_M',
'Y_SHIFT_M', 'ABS_SHIFT', 'ANGLE', 'SSIM_BEFORE', 'SSIM_AFTER',
'SSIM_IMPROVED', 'LAST_ERR'])
GDF = GDF.merge(records, on='POINT_ID', how="inner")
GDF = GDF.fillna(int(self.outFillVal))
......@@ -284,6 +305,13 @@ class Geom_Quality_Grid(object):
if getattr(GDF,'empty'): # GDF.empty returns AttributeError
return []
else:
orig_len_GDF = len(GDF)
if self.tieP_filter_level>0:
# SSIM filtering
GDF = GDF[GDF.SSIM_IMPROVED==True].copy()
if not self.q:
print('SSIM filtering removed %s tie points.' %(orig_len_GDF-len(GDF)))
# calculate GCPs
GDF['X_UTM_new'] = GDF.X_UTM + GDF.X_SHIFT_M
GDF['Y_UTM_new'] = GDF.Y_UTM + GDF.Y_SHIFT_M
......@@ -292,7 +320,7 @@ class Geom_Quality_Grid(object):
self.GCPList = GDF.GCP.tolist()
if not self.q:
print('Found %s valid GCPs.' %len(self.GCPList))
print('Found %s valid tie points.' %len(self.GCPList))
return self.GCPList
......
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