Commit 868d070e authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

implemented support for cloud masks / bad data masks

components.CoReg.imParams:
- __init__():  added attribute mask_baddata
- added add_mask_bad_data()

components.CoReg.COREG:
- __init__(): implemented new keywords 'mask_baddata_ref' and 'mask_baddata_tgt'
- _get_opt_winpos_winsize(): now checks bad data mask if window position is within bad data area

components.CoReg_local.COREG_LOCAL:
- __init__(): implemented new keywords 'mask_baddata_ref' and 'mask_baddata_tgt'

components.Geom_Quality_Grid.Geom_Quality_Grid:
- get_CoRegPoints_table(): added vectorized filtering of points within bad data area
- _flag_outliers(): added docstring, some minor revisions

coreg_cmd:
- renamed parameters path_im0 and path_im1 to path_ref and path_tgt
- added parameters '-mask_ref' and '-mask_tgt' to global an local argument parser

- updated __version__
parent f9bfdbb2
......@@ -9,7 +9,7 @@ from .components import utilities
from .components import geometry
__author__ = 'Daniel Scheffler'
__version__= '2016-11-12_01'
__version__= '2016-11-15_01'
__all__=['COREG',
'COREG_LOCAL',
......
......@@ -33,7 +33,7 @@ from py_tools_ds.ptds.geo.vector.topology import get_overlap_polygon, get_small
from py_tools_ds.ptds.geo.projection import prj_equal, get_proj4info
from py_tools_ds.ptds.geo.vector.geometry import boxObj, round_shapelyPoly_coords
from py_tools_ds.ptds.geo.coord_grid import move_shapelyPoly_to_image_grid
from py_tools_ds.ptds.geo.coord_trafo import pixelToMapYX, reproject_shapelyGeometry
from py_tools_ds.ptds.geo.coord_trafo import pixelToMapYX, reproject_shapelyGeometry, mapXY2imXY
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
......@@ -112,6 +112,33 @@ class imParamObj(object):
if not CoReg_params['q']:
print('Bounding box of calculated footprint for %s:\n\t%s' % (self.imName, self.poly.bounds))
# add bad data mask
given_mask = CoReg_params['mask_baddata_%s' % ('ref' if imID == 'ref' else 'tgt')]
self.mask_baddata = None
if given_mask:
self.add_mask_bad_data(given_mask)
def add_mask_bad_data(self, path_or_geoArr):
"""Adds a bad data mask. This method is separated from __init__() in order to allow explicit adding of the mask.
:param path_or_geoArr:
"""
geoArr_mask = path_or_geoArr if isinstance(path_or_geoArr, GeoArray) else GeoArray(path_or_geoArr)
assert geoArr_mask.bands == 1, \
'Expected one single band as bad data mask for %s. Got %s bands.' % (self.imName, geoArr_mask.bands)
pixelVals_in_mask = sorted(list(np.unique(geoArr_mask[:])))
assert len(pixelVals_in_mask) <= 2, 'Bad data mask must have only two pixel values (boolean) - 0 and 1 or ' \
'False and True! The given mask for %s contains the values %s.' % (
self.imName, pixelVals_in_mask)
assert pixelVals_in_mask in [[0, 1], [False, True]], 'Found unsupported pixel values in the given bad data ' \
'mask for %s: %s Only the values True, False, 0 and 1 ' \
'are supported. ' % (self.imName, pixelVals_in_mask)
self.mask_baddata = GeoArray(geoArr_mask[:].astype(np.bool), geoArr_mask.gt, geoArr_mask.prj)
class COREG(object):
......@@ -121,8 +148,9 @@ class COREG(object):
wp=(None,None), ws=(512, 512), max_iter=5, max_shift=5, align_grids=False, 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,
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):
nodata=(None,None), calc_corners=True, binary_ws=True, mask_baddata_ref=None, mask_baddata_tgt=None,
multiproc=True, force_quadratic_win=True, progress=True, v=False, path_verbose_out=None, q=False,
ignore_errors=False):
"""Detects and corrects global X/Y shifts between a target and refernce image. Geometric shifts are calculated
at a specific (adjustable) image position. Correction performs a global shifting in X- or Y direction.
......@@ -174,8 +202,20 @@ class COREG(object):
: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
(default: 1; deactivated if '-cor0' and '-cor1' are given
:param multiproc(bool): enable multiprocessing (default: 1)
:param binary_ws(bool): use binary X/Y dimensions for the matching window (default: 1)
:param mask_baddata_ref(str, GeoArray): path to a 2D boolean mask file (or an instance of GeoArray) for the
reference image where all bad data pixels (e.g. clouds) are marked with
True and the remaining pixels with False. Must have the same geographic
extent and projection like 'im_ref'. The mask is used to check if the
chosen matching window position is valid in the sense of useful data.
Otherwise this window position is rejected.
:param mask_baddata_tgt(str, GeoArray): path to a 2D boolean mask file (or an instance of GeoArray) for the
image to be shifted where all bad data pixels (e.g. clouds) are marked
with True and the remaining pixels with False. Must have the same
geographic extent and projection like 'im_ref'. The mask is used to
check if the chosen matching window position is valid in the sense of
useful data. Otherwise this window position is rejected.
:param multiproc(bool): enable multiprocessing (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: False)
......@@ -405,16 +445,20 @@ class COREG(object):
'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:
# if im.GeoArray[int(imY), int(imX), im.band4match]!=im.nodata:
# self.success = False
#else:
# assert im.GeoArray[int(imY), int(imX), im.band4match]!=im.nodata,\
# 'The provided window position is within the nodata area of the %s. Check the coordinates.' %im.imName
# check if window position is within bad data area if a respective mask has been provided
for im in [self.ref, self.shift]:
if im.mask_baddata is not None:
imX, imY = mapXY2imXY(wp, im.mask_baddata.gt)
if not im.mask_baddata[imY, imX]:
self.tracked_errors.append(
RuntimeError('According to the provided bad data mask for the %s the chosen window position '
'%s / %s is within a bad data area. Using this window position for coregistration '
'is not reasonable. Please provide a better window position!' %(im.imName, wp[0], wp[1])))
self.success = False
if not self.ignErr:
raise self.tracked_errors[-1]
self.win_pos_XY = wp
self.win_size_XY = (int(self.win_size_XY[0]), int(self.win_size_XY[1])) if self.win_size_XY else (512,512)
......
......@@ -29,7 +29,8 @@ class COREG_LOCAL(object):
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):
binary_ws=True, mask_baddata_ref=None, mask_baddata_tgt=None, 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
......@@ -93,6 +94,18 @@ class COREG_LOCAL(object):
matching window position within the actual image overlap
(default: True; deactivated if 'data_corners_im0' and 'data_corners_im1' are given
:param binary_ws(bool): use binary X/Y dimensions for the matching window (default: True)
:param mask_baddata_ref(str, GeoArray): path to a 2D boolean mask file (or an instance of GeoArray) for the
reference image where all bad data pixels (e.g. clouds) are marked with
True and the remaining pixels with False. Must have the same geographic
extent and projection like 'im_ref'. The mask is used to check if the
chosen matching window position is valid in the sense of useful data.
Otherwise this window position is rejected.
:param mask_baddata_tgt(str, GeoArray): path to a 2D boolean mask file (or an instance of GeoArray) for the
image to be shifted where all bad data pixels (e.g. clouds) are marked
with True and the remaining pixels with False. Must have the same
geographic extent and projection like 'im_ref'. The mask is used to
check if the chosen matching window position is valid in the sense of
useful data. Otherwise this window position is rejected.
:param CPUs(int): number of CPUs to use during calculation of geometric quality grid
(default: None, which means 'all CPUs available')
:param progress(bool): show progress bars (default: True)
......@@ -169,6 +182,11 @@ class COREG_LOCAL(object):
q = q,
ignore_errors = self.ignErr)
# add bad data mask
# (mask is not added during initialization of COREG object in order to avoid bad data area errors there)
if mask_baddata_ref: self.COREG_obj.ref.add_mask_bad_data(mask_baddata_ref)
if mask_baddata_tgt: self.COREG_obj.shift.add_mask_bad_data(mask_baddata_tgt)
self._quality_grid = None # set by self.quality_grid
self._CoRegPoints_table = None # set by self.CoRegPoints_table
self._coreg_info = None # set by self.coreg_info
......
......@@ -16,7 +16,8 @@ import numpy as np
from geopandas import GeoDataFrame
from pykrige.ok import OrdinaryKriging
from shapely.geometry import Point
from skimage.measure import points_in_poly
from skimage.measure import points_in_poly, ransac
from skimage.transform import AffineTransform, PolynomialTransform
# internal modules
from .CoReg import COREG
......@@ -167,7 +168,7 @@ class Geom_Quality_Grid(object):
return [pointID]+CR_res
def get_CoRegPoints_table(self, exclude_outliers=1):
def get_CoRegPoints_table(self, exclude_outliers=True):
assert self.XY_points is not None and self.XY_mapPoints is not None
#ref_ds,tgt_ds = gdal.Open(self.path_imref),gdal.Open(self.path_im2shift)
......@@ -182,7 +183,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]) # FIXME replace that with coord image reprojection
XYarr2PointGeom = np.vectorize(lambda X,Y: Point(X,Y), otypes=[Point])
geomPoints = np.array(XYarr2PointGeom(self.XY_mapPoints[:,0],self.XY_mapPoints[:,1]))
if isProjectedOrGeographic(self.COREG_obj.shift.prj)=='geographic':
......@@ -206,9 +207,20 @@ class Geom_Quality_Grid(object):
inliers = points_in_poly(self.XY_mapPoints,
np.swapaxes(np.array(self.COREG_obj.overlap_poly.exterior.coords.xy), 0, 1))
GDF = GDF[inliers]
#GDF = GDF[GDF['geometry'].within(self.COREG_obj.overlap_poly.simplify(tolerance=15))] # slow
#GDF = GDF[GDF['geometry'].within(self.COREG_obj.overlap_poly.simplify(tolerance=15))] # works but much slower
assert not GDF.empty, 'No coregistration point could be placed within the overlap area. Check yout input data!' # FIXME track that
assert not GDF.empty, 'No coregistration point could be placed within the overlap area. Check your input data!' # FIXME track that
# exclude all point where bad data mask is True (e.g. points on clouds etc.)
orig_len_GDF = len(GDF)
mapXY = np.array(GDF.loc[:,['X_UTM','Y_UTM']])
GDF['REF_BADDATA'] = self.COREG_obj.ref .mask_baddata.read_pointData(mapXY) \
if self.COREG_obj.ref .mask_baddata is not None else False
GDF['TGT_BADDATA'] = self.COREG_obj.shift.mask_baddata.read_pointData(mapXY)\
if self.COREG_obj.shift.mask_baddata is not None else False
GDF = GDF[(GDF['REF_BADDATA']==False) & (GDF['TGT_BADDATA']==False)]
print('According to the provided bad data mask(s) %s points of initially %s have been excluded.'
%(orig_len_GDF-len(GDF), orig_len_GDF))
#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 \
......@@ -288,25 +300,30 @@ class Geom_Quality_Grid(object):
self.CoRegPoints_table = GDF
if self.tieP_filter_level>1:
self._flag_outliers(algorithm = 'RANSAC')
self._flag_outliers(algorithm = 'RANSAC', inplace=True)
return self.CoRegPoints_table
def _flag_outliers(self, algorithm='RANSAC'):
def _flag_outliers(self, algorithm='RANSAC', inplace=False):
"""Detects geometric outliers within CoRegPoints_table GeoDataFrame and adds a new boolean column 'OUTLIER'.
:param algorithm: <str> the algorithm to be used for outlier detection. available choices: 'RANSAC'
:param inplace: <bool> whether to overwrite CoRegPoints_table directly (True) or not (False). If False, the
resulting GeoDataFrame is returned.
"""
if not algorithm in ['RANSAC']: raise ValueError
GDF = self.CoRegPoints_table
GDF = GDF[GDF.ABS_SHIFT!=self.outFillVal]
#GDF = GDF[(GDF.ABS_SHIFT!=self.outFillVal) &(GDF.SSIM_IMPROVED==True)]
GDF = self.CoRegPoints_table.copy()
# exclude all records where SSIM decreased or no match has been found
GDF = GDF[(GDF.ABS_SHIFT!=self.outFillVal) &(GDF.SSIM_IMPROVED==True)]
src = np.array(GDF[['X_UTM', 'Y_UTM']])
xyShift = np.array(GDF[['X_SHIFT_M', 'Y_SHIFT_M']])
dst = src + xyShift
if algorithm=='RANSAC':
from skimage.measure import ransac
from skimage.transform import AffineTransform, PolynomialTransform
# robustly estimate affine transform model with RANSAC
#model_robust, inliers = ransac((src, dst), PolynomialTransform, min_samples=3,
model_robust, inliers = ransac((src, dst), AffineTransform,
......@@ -320,7 +337,13 @@ class Geom_Quality_Grid(object):
records = GDF[['POINT_ID']].copy()
records['OUTLIER'] = outliers
self.CoRegPoints_table = self.CoRegPoints_table.merge(records, on='POINT_ID', how="inner")
GDF = self.CoRegPoints_table.merge(records, on='POINT_ID', how="outer")
GDF['OUTLIER'].fillna(int(self.outFillVal), inplace=True)
if inplace:
self.CoRegPoints_table = GDF
return GDF
def dump_CoRegPoints_table(self, path_out=None):
......
......@@ -39,8 +39,8 @@ except ImportError:
#sub-command functions
def run_global_coreg(args):
COREG_obj = COREG(args.path_im0,
args.path_im1,
COREG_obj = COREG(args.path_ref,
args.path_tgt,
path_out = args.path_out,
fmt_out = args.fmt_out,
r_b4match = args.br,
......@@ -58,6 +58,8 @@ def run_global_coreg(args):
calc_corners = args.calc_cor,
multiproc = args.mp,
binary_ws = args.bin_ws,
mask_baddata_ref = args.mask_ref,
mask_baddata_tgt = args.mask_tgt,
v = args.v,
path_verbose_out = args.vo,
q = args.q,
......@@ -67,8 +69,8 @@ def run_global_coreg(args):
#sub-command functions
def run_local_coreg(args):
CRL = COREG_LOCAL(args.path_im0,
args.path_im1,
CRL = COREG_LOCAL(args.path_ref,
args.path_tgt,
path_out = args.path_out,
fmt_out = args.fmt_out,
grid_res = args.grid_res,
......@@ -84,6 +86,8 @@ def run_local_coreg(args):
data_corners_tgt = args.cor1,
nodata = args.nodata,
calc_corners = args.calc_cor,
mask_baddata_ref = args.mask_ref,
mask_baddata_tgt = args.mask_tgt,
CPUs = None if args.mp else 1,
binary_ws = args.bin_ws,
progress = args.progress,
......@@ -150,9 +154,9 @@ if __name__ == '__main__':
"use '>>> python /path/to/CoReg_Sat/coreg_cmd.py global -h' for documentation and usage hints")
gloArg = parse_coreg_global.add_argument
gloArg('path_im0', type=str, help='source path of reference image (any GDAL compatible image format is supported)')
gloArg('path_ref', type=str, help='source path of reference image (any GDAL compatible image format is supported)')
gloArg('path_im1', type=str, help='source path of image to be shifted (any GDAL compatible image format is supported)')
gloArg('path_tgt', type=str, help='source path of image to be shifted (any GDAL compatible image format is supported)')
gloArg('-o', nargs='?', dest='path_out', type=str, default='auto',
help="target path of the coregistered image (default: /dir/of/im1/<im1>__shifted_to__<im0>.bsq)")
......@@ -204,14 +208,26 @@ if __name__ == '__main__':
help="calculate true positions of the dataset corners in order to get a useful matching window position "
"within the actual image overlap (default: 1; deactivated if '-cor0' and '-cor1' are given")
gloArg('-mp', nargs='?', type=int, help='enable multiprocessing (default: 1)', default=1, choices=[0, 1])
gloArg('-bin_ws', nargs='?', type=int,
help='use binary X/Y dimensions for the matching window (default: 1)', default=1, choices=[0, 1])
gloArg('-quadratic_win', nargs='?', type=int,
help='force a quadratic matching window (default: 1)', default=1, choices=[0, 1])
gloArg('-mask_ref', nargs='?', type=str, default=None, metavar='file path',
help="path to a 2D boolean mask file for the reference image where all bad data pixels (e.g. clouds) are "
"marked with True or 1 and the remaining pixels with False or 0. Must have the same geographic extent "
"and projection like the refernce image. The mask is used to check if the chosen matching window "
"position is valid in the sense of useful data. Otherwise this window position is rejected.")
gloArg('-mask_tgt', nargs='?', type=str, default=None, metavar='file path',
help="path to a 2D boolean mask file for the image to be shifted where all bad data pixels (e.g. clouds) are "
"marked with True or 1 and the remaining pixels with False or 0. Must have the same geographic extent "
"and projection like the the image to be shifted. The mask is used to check if the chosen matching "
"window position is valid in the sense of useful data. Otherwise this window position is rejected.")
gloArg('-mp', nargs='?', type=int, help='enable multiprocessing (default: 1)', default=1, choices=[0, 1])
gloArg('-v', nargs='?', type=int, help='verbose mode (default: 0)', default=0, choices=[0, 1])
gloArg('-vo', nargs='?', type=str, choices=[0, 1], help='an optional output directory for outputs of verbose mode'
......@@ -237,9 +253,9 @@ if __name__ == '__main__':
"use '>>> python /path/to/CoReg_Sat/coreg_cmd.py local -h' for documentation and usage hints")
locArg = parse_coreg_local.add_argument
locArg('path_im0', type=str, help='source path of reference image (any GDAL compatible image format is supported)')
locArg('path_ref', type=str, help='source path of reference image (any GDAL compatible image format is supported)')
locArg('path_im1', type=str, help='source path of image to be shifted (any GDAL compatible image format is supported)')
locArg('path_tgt', type=str, help='source path of image to be shifted (any GDAL compatible image format is supported)')
locArg('grid_res', type=int, help='quality grid resolution in pixels of the target image')
......@@ -287,6 +303,18 @@ if __name__ == '__main__':
locArg('-quadratic_win', nargs='?', type=int,
help='force a quadratic matching window (default: 1)', default=1, choices=[0, 1])
locArg('-mask_ref', nargs='?', type=str, default=None, metavar='file path',
help="path to a 2D boolean mask file for the reference image where all bad data pixels (e.g. clouds) are "
"marked with True or 1 and the remaining pixels with False or 0. Must have the same geographic extent "
"and projection like the refernce image. The mask is used to check if the chosen matching window "
"position is valid in the sense of useful data. Otherwise this window position is rejected.")
locArg('-mask_tgt', nargs='?', type=str, default=None, metavar='file path',
help="path to a 2D boolean mask file for the image to be shifted where all bad data pixels (e.g. clouds) are "
"marked with True or 1 and the remaining pixels with False or 0. Must have the same geographic extent "
"and projection like the the image to be shifted. The mask is used to check if the chosen matching "
"window position is valid in the sense of useful data. Otherwise this window position is rejected.")
locArg('-progress', nargs='?', type=int, help='show progress bars (default: 1)', default=1, choices=[0, 1])
locArg('-v', nargs='?', type=int, help='verbose mode (default: 0)', default=0, choices=[0, 1])
......
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