Commit 57487d9d authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'dev'

parents d276bbf8 9736e3b6
......@@ -19,7 +19,7 @@ try:
import pyfftw
except ImportError:
pyfftw = None
from shapely.geometry import Point
from shapely.geometry import Point, Polygon
# internal modules
from .DeShifter import DESHIFTER, _dict_rspAlg_rsp_Int
......@@ -87,32 +87,29 @@ class imParamObj(object):
else:
self.nodata = self.GeoArray.nodata
# set corner coords
given_corner_coord = CoReg_params['data_corners_%s' % ('im0' if imID == 'ref' else 'im1')]
if given_corner_coord is None:
if CoReg_params['calc_corners']:
if not CoReg_params['q']:
print('Calculating actual data corner coordinates for %s...' % self.imName)
self.corner_coord = GEO.get_true_corner_mapXY(self.GeoArray, self.band4match, self.nodata,
CoReg_params['multiproc'], v=self.v, q=self.q)
else:
self.corner_coord = get_corner_coordinates(gt=self.GeoArray.geotransform,
cols=self.cols,rows=self.rows)
# set footprint_poly
given_footprint_poly = CoReg_params['footprint_poly_%s' % ('ref' if imID == 'ref' else 'tgt')]
given_corner_coord = CoReg_params['data_corners_%s' % ('ref' if imID == 'ref' else 'tgt')]
if given_footprint_poly:
self.GeoArray.footprint_poly = given_footprint_poly
elif given_corner_coord is not None:
self.GeoArray.footprint_poly = Polygon(given_corner_coord)
elif not CoReg_params['calc_corners']:
# use the image extent
self.GeoArray.footprint_poly = Polygon(get_corner_coordinates(gt=self.GeoArray.geotransform,
cols=self.cols,rows=self.rows))
else:
self.corner_coord = given_corner_coord
# set footprint polygon
#self.poly = get_footprint_polygon(self.corner_coord, fix_invalid=True) # this is the old algorithm
self.GeoArray.calc_mask_nodata(fromBand=self.band4match) # this avoids that all bands have to be read
# footprint_poly is calculated automatically by GeoArray
if not CoReg_params['q']:
print('Calculating actual data corner coordinates for %s...' % self.imName)
self.GeoArray.calc_mask_nodata(fromBand=self.band4match) # this avoids that all bands have to be read
self.poly = self.GeoArray.footprint_poly
self.poly = self.GeoArray.footprint_poly # returns a shapely geometry
for XY in self.corner_coord:
assert self.GeoArray.box.mapPoly.contains(Point(XY)) or self.GeoArray.box.mapPoly.touches(Point(XY)), \
"The corner position '%s' is outside of the %s." % (XY, self.imName)
if not CoReg_params['q']:
print('Bounding box of calculated footprint for %s:\n\t%s' % (self.imName, self.poly.bounds))
if not CoReg_params['q']: print('Corner coordinates of %s:\n\t%s' % (self.imName, self.corner_coord))
class COREG(object):
......@@ -120,7 +117,8 @@ class COREG(object):
def __init__(self, im_ref, im_tgt, path_out=None, fmt_out='ENVI', r_b4match=1, s_b4match=1, wp=(None,None),
ws=(512, 512), max_iter=5, max_shift=5, align_grids=False, match_gsd=False, out_gsd=None,
resamp_alg_deshift='cubic', resamp_alg_calc='cubic', data_corners_im0=None, data_corners_im1=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,
nodata=(None,None), calc_corners=True, multiproc=True, binary_ws=True, force_quadratic_win=True,
progress=True, v=False, path_verbose_out=None, q=False, ignore_errors=False):
......@@ -156,8 +154,10 @@ class COREG(object):
(valid algorithms: nearest, bilinear, cubic, cubic_spline, lanczos, average, mode,
max, min, med, q1, q3)
default: cubic (highly recommended)
:param data_corners_im0(list): map coordinates of data corners within reference image
:param data_corners_im1(list): map coordinates of data corners within image to be shifted
:param footprint_poly_ref(str):
:param footprint_poly_tgt(str):
:param data_corners_ref(list): map coordinates of data corners within reference image
:param data_corners_tgt(list): map coordinates of data corners within image to be shifted
:param nodata(tuple): no data values for reference image and image to be shifted
:param calc_corners(bool): calculate true positions of the dataset corners in order to get a useful
matching window position within the actual image overlap
......@@ -166,11 +166,11 @@ class COREG(object):
:param binary_ws(bool): use binary X/Y dimensions for the matching window (default: 1)
:param force_quadratic_win(bool): force a quadratic matching window (default: 1)
:param progress(bool): show progress bars (default: True)
:param v(bool): verbose mode (default: 0)
:param v(bool): verbose mode (default: False)
:param path_verbose_out(str): an optional output directory for intermediate results
(if not given, no intermediate results are written to disk)
:param q(bool): quiet mode (default: 0)
:param ignore_errors(bool): Useful for batch processing. (default: 0)
:param q(bool): quiet mode (default: False)
:param ignore_errors(bool): Useful for batch processing. (default: False)
In case of error COREG.success == False and COREG.x_shift_px/COREG.y_shift_px
is None
"""
......@@ -179,10 +179,10 @@ class COREG(object):
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.'
if data_corners_im0 and not isinstance(data_corners_im0[0],list): # group if not [[x,y],[x,y]..] but [x,y,x,y,]
data_corners_im0 = [data_corners_im0[i:i+2] for i in range(0, len(data_corners_im0), 2)]
if data_corners_im1 and not isinstance(data_corners_im1[0],list): # group if not [[x,y],[x,y]..]
data_corners_im1 = [data_corners_im1[i:i+2] for i in range(0, len(data_corners_im1), 2)]
if data_corners_ref and not isinstance(data_corners_ref[0], list): # group if not [[x,y],[x,y]..] but [x,y,x,y,]
data_corners_ref = [data_corners_ref[i:i + 2] for i in range(0, len(data_corners_ref), 2)]
if data_corners_tgt and not isinstance(data_corners_tgt[0], list): # group if not [[x,y],[x,y]..]
data_corners_tgt = [data_corners_tgt[i:i + 2] for i in range(0, len(data_corners_tgt), 2)]
if nodata: assert isinstance(nodata, tuple) and len(nodata) == 2, "'nodata' must be a tuple with two values." \
"Got %s with length %s." %(type(nodata),len(nodata))
for rspAlg in [resamp_alg_deshift, resamp_alg_calc]:
......@@ -230,7 +230,7 @@ class COREG(object):
self.updated_map_info = None # set by self.get_updated_map_info()
self.tracked_errors = [] # expanded each time an error occurs
self.success = False # default
self.success = None # default
self.deshift_results = None # set by self.correct_shifts()
gdal.AllRegister()
......@@ -257,11 +257,11 @@ class COREG(object):
if not self.q: print('Matching window position (X,Y): %s/%s' % (self.win_pos_XY[0], self.win_pos_XY[1]))
self._get_clip_window_properties()
if self.v and self.path_verbose_out and self.matchWin.mapPoly:
if self.v and self.path_verbose_out and self.matchWin.mapPoly and self.success is not False:
IO.write_shp(os.path.join(self.path_verbose_out, 'poly_matchWin.shp'),
self.matchWin.mapPoly, self.matchWin.prj)
self.success = None if self.matchWin.boxMapYX else False
self.success = False if self.success is False or not self.matchWin.boxMapYX else None
self._coreg_info = None # private attribute to be filled by self.coreg_info property
......@@ -438,23 +438,36 @@ class COREG(object):
# -> match Fenster verkleinern und neues anderes Fenster berechnen
xLarger, yLarger = otherWin.is_larger_DimXY(overlapWin.boundsIm)
matchWin.buffer_imXY(-1 if xLarger else 0, -1 if yLarger else 0)
previous_area = otherWin.mapPoly.area
otherWin.boxImYX = get_smallest_boxImYX_that_contains_boxMapYX(matchWin.boxMapYX,otherWin.gt)
# check results
assert matchWin.mapPoly.within(otherWin.mapPoly)
assert otherWin.mapPoly.within(overlapWin.mapPoly)
self.imfft_gsd = self.ref.xgsd if self.grid2use =='ref' else self.shift.xgsd
self.ref.win,self.shift.win = (matchWin,otherWin) if self.grid2use =='ref' else (otherWin,matchWin)
self.matchWin,self.otherWin = matchWin, otherWin
self.ref. win.size_YX = tuple([int(i) for i in self.ref. win.imDimsYX])
self.shift.win.size_YX = tuple([int(i) for i in self.shift.win.imDimsYX])
match_win_size_XY = tuple(reversed([int(i) for i in matchWin.imDimsYX]))
if not self.q and match_win_size_XY != self.win_size_XY:
print('Target window size %s not possible due to too small overlap area or window position too close '
'to an image edge. New matching window size: %s.' %(self.win_size_XY,match_win_size_XY))
#IO.write_shp('/misc/hy5/scheffler/Temp/matchMapPoly.shp', matchWin.mapPoly,matchWin.prj)
#IO.write_shp('/misc/hy5/scheffler/Temp/otherMapPoly.shp', otherWin.mapPoly,otherWin.prj)
if previous_area == otherWin.mapPoly.area:
self.tracked_errors.append(
RuntimeError('Matching window in target image is larger than overlap area but further shrinking '
'the matching window is not possible. Check if the footprints of the input data have '
'been computed correctly. '))
if not self.ignErr:
raise self.tracked_errors[-1]
break # break out of while loop in order to avoid that code gets stuck here
if self.tracked_errors:
self.success = False
else:
# check results
assert matchWin.mapPoly.within(otherWin.mapPoly)
assert otherWin.mapPoly.within(overlapWin.mapPoly)
self.imfft_gsd = self.ref.xgsd if self.grid2use =='ref' else self.shift.xgsd
self.ref.win,self.shift.win = (matchWin,otherWin) if self.grid2use =='ref' else (otherWin,matchWin)
self.matchWin,self.otherWin = matchWin, otherWin
self.ref. win.size_YX = tuple([int(i) for i in self.ref. win.imDimsYX])
self.shift.win.size_YX = tuple([int(i) for i in self.shift.win.imDimsYX])
match_win_size_XY = tuple(reversed([int(i) for i in matchWin.imDimsYX]))
if not self.q and match_win_size_XY != self.win_size_XY:
print('Target window size %s not possible due to too small overlap area or window position too close '
'to an image edge. New matching window size: %s.' %(self.win_size_XY,match_win_size_XY))
#IO.write_shp('/misc/hy5/scheffler/Temp/matchMapPoly.shp', matchWin.mapPoly,matchWin.prj)
#IO.write_shp('/misc/hy5/scheffler/Temp/otherMapPoly.shp', otherWin.mapPoly,otherWin.prj)
def _get_image_windows_to_match(self):
......
......@@ -25,9 +25,9 @@ class COREG_LOCAL(object):
"""See help(COREG_LOCAL) for documentation!"""
def __init__(self, im_ref, im_tgt, grid_res, window_size=(256,256), path_out=None, fmt_out='ENVI', projectDir=None,
r_b4match=1, s_b4match=1, max_iter=5, max_shift=5, data_corners_im0=None,
data_corners_im1=None, outFillVal=-9999, nodata=(None,None), calc_corners=True, binary_ws=True,
CPUs=None, progress=True, v=False, q=False):
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):
"""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
......@@ -48,8 +48,14 @@ 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 data_corners_im0(list): map coordinates of data corners within reference image
:param data_corners_im1(list): map coordinates of data corners within image to be shifted
: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))'
:param footprint_poly_tgt(str): footprint polygon of the image to be shifted (WKT string or shapely.geometry.Polygon)
e.g. 'POLYGON ((299999 6000000, 299999 5890200, 409799 5890200, 409799 6000000,
299999 6000000))'
:param data_corners_ref(list): map coordinates of data corners within reference image
:param data_corners_tgt(list): map coordinates of data corners within image to be shifted
: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 nodata(tuple): no data values for reference image and image to be shifted
......@@ -62,6 +68,7 @@ class COREG_LOCAL(object):
:param progress(bool): show progress bars (default: True)
:param v(bool): verbose mode (default: False)
:param q(bool): quiet mode (default: False)
:param ignore_errors(bool): Useful for batch processing. (default: False)
"""
# TODO outgsd, matchgsd
# assertions
......@@ -95,6 +102,7 @@ class COREG_LOCAL(object):
self.v = v
self.q = q if not v else False # overridden by v
self.progress = progress if not q else False # overridden by v
self.ignErr = ignore_errors
COREG.__dict__['_set_outpathes'](self, self.imref, self.im2shift)
# make sure that the output directory of coregistered image is the project directory if a project directory is given
......@@ -103,21 +111,23 @@ class COREG_LOCAL(object):
gdal.AllRegister()
self.COREG_obj = COREG(self.imref, self.im2shift,
ws = window_size,
data_corners_im0 = data_corners_im0,
data_corners_im1 = data_corners_im1,
calc_corners = calc_corners,
r_b4match = r_b4match,
s_b4match = s_b4match,
max_iter = max_iter,
max_shift = max_shift,
nodata = nodata,
multiproc = self.CPUs is None or self.CPUs > 1,
binary_ws = self.bin_ws,
progress = self.progress,
v = v,
q = q,
ignore_errors = True)
ws = window_size,
footprint_poly_ref = footprint_poly_ref,
footprint_poly_tgt = footprint_poly_tgt,
data_corners_ref = data_corners_ref,
data_corners_tgt = data_corners_tgt,
calc_corners = calc_corners,
r_b4match = r_b4match,
s_b4match = s_b4match,
max_iter = max_iter,
max_shift = max_shift,
nodata = nodata,
multiproc = self.CPUs is None or self.CPUs > 1,
binary_ws = self.bin_ws,
progress = self.progress,
v = v,
q = q,
ignore_errors = self.ignErr)
self._quality_grid = None # set by self.quality_grid
self._CoRegPoints_table = None # set by self.CoRegPoints_table
......
......@@ -180,6 +180,8 @@ class Geom_Quality_Grid(object):
if exclude_outliers:
GDF = GDF[GDF['geometry'].within(self.COREG_obj.overlap_poly)]
assert not GDF.empty, 'No coregistration point could be placed within the overlap area. Check yout input data!' # FIXME track that
#not_within_nodata = \
# lambda r: np.array(self.ref[r.Y_IM,r.X_IM,self.COREG_obj.ref.band4match]!=self.COREG_obj.ref.nodata and \
# self.shift[r.Y_IM,r.X_IM, self.COREG_obj.shift.band4match] != self.COREG_obj.shift.nodata)[0,0]
......@@ -198,20 +200,20 @@ class Geom_Quality_Grid(object):
# get all variations of kwargs for coregistration
get_coreg_kwargs = lambda pID, wp: {
'pointID' : pID,
'wp' : wp,
'ws' : self.COREG_obj.win_size_XY,
'data_corners_im0': self.COREG_obj.ref.corner_coord,
'data_corners_im1': self.COREG_obj.shift.corner_coord,
'r_b4match' : self.COREG_obj.ref.band4match+1, # band4match is internally saved as index, starting from 0
's_b4match' : self.COREG_obj.shift.band4match+1, # band4match is internally saved as index, starting from 0
'max_iter' : self.COREG_obj.max_iter,
'max_shift' : self.COREG_obj.max_shift,
'nodata' : (self.COREG_obj.ref.nodata, self.COREG_obj.shift.nodata),
'binary_ws' : self.COREG_obj.bin_ws,
'v' : False, # otherwise this would lead to massive console output
'q' : True, # otherwise this would lead to massive console output
'ignore_errors' : True
'pointID' : pID,
'wp' : wp,
'ws' : self.COREG_obj.win_size_XY,
'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
's_b4match' : self.COREG_obj.shift.band4match+1, # band4match is internally saved as index, starting from 0
'max_iter' : self.COREG_obj.max_iter,
'max_shift' : self.COREG_obj.max_shift,
'nodata' : (self.COREG_obj.ref.nodata, self.COREG_obj.shift.nodata),
'binary_ws' : self.COREG_obj.bin_ws,
'v' : False, # otherwise this would lead to massive console output
'q' : True, # otherwise this would lead to massive console output
'ignore_errors' : True
}
list_coreg_kwargs = (get_coreg_kwargs(i, self.XY_mapPoints[i]) for i in GDF.index) # generator
......
......@@ -46,6 +46,7 @@ def get_true_corner_mapXY(fPath_or_geoarray, band=0, noDataVal=None, mp=1, v=0,
:param q:
:return:
"""
warnings.warn('This function is not in use anymore. Use it on your own risk!', DeprecationWarning)
geoArr = GeoArray(fPath_or_geoarray) if not isinstance(fPath_or_geoarray,GeoArray) else fPath_or_geoarray
rows,cols = geoArr.shape[:2]
......
......@@ -51,8 +51,8 @@ def run_global_coreg(args):
align_grids = args.align_grids,
match_gsd = args.match_gsd,
out_gsd = args.out_gsd,
data_corners_im0 = args.cor0,
data_corners_im1 = args.cor1,
data_corners_ref = args.cor0,
data_corners_tgt = args.cor1,
nodata = args.nodata,
calc_corners = args.calc_cor,
multiproc = args.mp,
......@@ -79,8 +79,8 @@ def run_local_coreg(args):
#align_grids = args.align_grids,
#match_gsd = args.match_gsd,
#out_gsd = args.out_gsd,
data_corners_im0 = args.cor0,
data_corners_im1 = args.cor1,
data_corners_ref = args.cor0,
data_corners_tgt = args.cor1,
nodata = args.nodata,
calc_corners = args.calc_cor,
CPUs = None if args.mp else 1,
......@@ -183,6 +183,8 @@ if __name__ == '__main__':
gloArg('-out_gsd', nargs=2, type=float, help='xgsd ygsd: set the output pixel size in map units'\
'(default: original pixel size of the image to be shifted)', metavar=('xgsd','ygsd'))
# TODO implement footprint_poly_ref, footprint_poly_tgt
gloArg('-cor0', nargs=8, type=float, help="map coordinates of data corners within reference image: ",
metavar=tuple("UL-X UL-Y UR-X UR-Y LR-X LR-Y LL-X LL-Y".split(' ')), default=None)
......@@ -256,6 +258,8 @@ if __name__ == '__main__':
locArg('-max_shift', nargs='?', type=int,
help="maximum shift distance in reference image pixel units (default: 5 px)", default=5)
# TODO implement footprint_poly_ref, footprint_poly_tgt
locArg('-cor0', nargs=8, type=float, help="map coordinates of data corners within reference image: ",
metavar=tuple("UL-X UL-Y UR-X UR-Y LR-X LR-Y LL-X LL-Y".split(' ')), default=None)
......
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