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


- imParamObj: simplified parameter settings and adjusted calls of GEO.find_noDataVal

- get_quality_grid(): bugfix for processing the wrong input array shape during coregistration
- added quiet mode for some print outputs

- added new version of find_noDataVal() that now also supports GeoArrays as input
- get_true_corner_mapXY(): added automatic reordering of corner coordinates in order to avoid shapely self intersection warnings
parent 240aea01
......@@ -39,6 +39,7 @@ from py_tools_ds.ptds.numeric.vector import find_nearest
class imParamObj(object):
def __init__(self, CoReg_params, imID):
assert imID in ['ref', 'shift']
self.imName = 'reference image' if imID == 'ref' else 'image to be shifted'
self.v = CoReg_params['v']
self.q = CoReg_params['q'] if not self.v else False
......@@ -53,12 +54,11 @@ class imParamObj(object):
# set params
self.prj = self.GeoArray.projection = self.GeoArray.geotransform
self.xgsd = abs([1])
self.ygsd = abs([5])
shape = self.GeoArray.shape
self.rows = shape[0]
self.cols = shape[1]
self.bands = shape[2] if len(shape)==3 else 1
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
# validate params
assert self.prj, 'The %s has no projection.' % self.imName
......@@ -72,14 +72,10 @@ class imParamObj(object):
if self.bands > 1 else '', self.bands)
# set nodata
if CoReg_params['nodata'][0] is not None and CoReg_params['nodata'][1] is not None:
if CoReg_params['nodata'][0 if imID == 'ref' else 1] is not None:
self.nodata = CoReg_params['nodata'][0 if imID == 'ref' else 1]
elif self.GeoArray.filePath:
self.nodata = GEO.find_noDataVal(self.GeoArray.filePath)
else: #FIXME
warnings.warn('Automatic nodata value detection for numpy array is currently not implemented. '
'Please provide a no data value. Otherwise it is set to None.')
self.nodata = None
self.nodata = GEO.find_noDataVal(self.GeoArray)
# set corner coords
given_corner_coord = CoReg_params['data_corners_%s' % ('im0' if imID == 'ref' else 'im1')]
......@@ -180,9 +180,9 @@ class Geom_Quality_Grid(object):
# declare global variables needed for self._get_spatial_shifts()
global global_shared_imref,global_shared_im2shift
global_shared_imref = \
GeoArray(self.imref[self.r_b4match-1], self.imref.geotransform, self.imref.projection)
GeoArray(self.imref[:,:,self.r_b4match-1], self.imref.geotransform, self.imref.projection)
global_shared_im2shift = \
GeoArray(self.im2shift[self.s_b4match-1], self.im2shift.geotransform,self.im2shift.projection)
GeoArray(self.im2shift[:,:,self.s_b4match-1], self.im2shift.geotransform,self.im2shift.projection)
# get all variations of kwargs for coregistration
get_coreg_kwargs = lambda pID, wp: {
......@@ -205,11 +205,13 @@ class Geom_Quality_Grid(object):
# run co-registration for whole grid
if not self.q:
print("Calculating geometric quality grid in mode 'multiprocessing'...")
with multiprocessing.Pool() as pool:
results =, list_coreg_kwargs)
if not self.q:
print("Calculating geometric quality grid in mode 'multiprocessing'...")
results = []
for i,coreg_kwargs in enumerate(list_coreg_kwargs):
......@@ -288,7 +290,8 @@ class Geom_Quality_Grid(object):
os.path.basename(self.im2shift.filePath))[0], os.path.splitext(os.path.basename(self.imref.filePath))[0]) # FIXME does not work for inmem GeoArrays
path_out = os.path.join(self.dir_out, 'CoRegPoints', fName_out)
if not os.path.exists(os.path.dirname(path_out)): os.makedirs(os.path.dirname(path_out))
print('Writing %s ...' %path_out)
if not self.q:
print('Writing %s ...' %path_out)
......@@ -65,6 +65,9 @@ def get_true_corner_mapXY(fPath_or_geoarray, bandNr=1, noDataVal=None, mp=1, v=0
" Using algorithm 'numpy' instead." %sys.exc_info()[1])
corner_coords_YX = calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=False, algorithm='numpy')
if len(corner_coords_YX)==4: # this avoids shapely self intersection
corner_coords_YX = list(np.array(corner_coords_YX)[[0,1,3,2]]) # UL, UR, LL, LR => UL, UR, LR, LL
# check if enough unique coordinates have been found
if not len(GeoDataFrame(corner_coords_YX).drop_duplicates().values)>=3:
if not q:
......@@ -102,21 +105,18 @@ def get_GeoArrayPosition_from_boxImYX(boxImYX):
return rS, rE, cS, cE
def find_noDataVal(path_im,bandNr=1,sz=3,is_vrt=0):
def find_noDataVal(pathIm_or_GeoArray,bandIdx=0,sz=3):
"""tries to derive no data value from homogenious corner pixels within 3x3 windows (by default)
:param path_im:
:param bandNr:
:param pathIm_or_GeoArray:
:param bandIdx:
:param sz: window size in which corner pixels are analysed
:param is_vrt:
gdal_ds = gdal.Open(path_im) if not is_vrt else gdal.OpenShared(path_im)
geoArr = pathIm_or_GeoArray if isinstance(pathIm_or_GeoArray, GeoArray) else GeoArray(pathIm_or_GeoArray)
get_mean_std = lambda corner_subset: {'mean':np.mean(corner_subset), 'std':np.std(corner_subset)}
rows,cols = gdal_ds.RasterYSize,gdal_ds.RasterXSize
UL = get_mean_std(gdal_ds.GetRasterBand(bandNr).ReadAsArray(0,0,sz,sz))
UR = get_mean_std(gdal_ds.GetRasterBand(bandNr).ReadAsArray(cols-sz,0,sz,sz))
LR = get_mean_std(gdal_ds.GetRasterBand(bandNr).ReadAsArray(cols-sz,rows-sz,sz,sz))
LL = get_mean_std(gdal_ds.GetRasterBand(bandNr).ReadAsArray(0,rows-sz,sz,sz))
gdal_ds = None
UL = get_mean_std(geoArr[0:sz,0:sz,bandIdx])
UR = get_mean_std(geoArr[0:sz,-sz:,bandIdx])
LR = get_mean_std(geoArr[-sz:,-sz:,bandIdx])
LL = get_mean_std(geoArr[-sz:,0:sz,bandIdx])
possVals = [i['mean'] for i in [UL,UR,LR,LL] if i['std']==0]
# possVals==[]: all corners are filled with data; np.std(possVals)==0: noDataVal clearly identified
return None if possVals==[] else possVals[0] if np.std(possVals)==0 else 'unclear'
\ No newline at end of file
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