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

bugfix for outdated assertion; new keyword to switch off progress bars, some further developments

components.CoReg:
- imParams: bugfix for outdated assertion
- COREG:
    - implemented keyword 'progress' for switching progress bars on or off
    - converted attribute 'coreg_info' to property

components.CoReg_local.COREG_LOCAL:
- implemented keyword 'progress' for switching progress bars on or off
- added property 'coreg_info'

components.DeShifter.DESHIFTER:
- implemented keyword 'progress' for switching progress bars on or off
- removed deprecated keywords 'warp_alg' and 'tempDir'
- correct_shifts(): outcommented algorithm for GDAL command line warping

components.Geom_Quality_Grid.Geom_Quality_Grid:
- implemented keyword 'progress' for switching progress bars on or off
- get_CoRegPoints_table():
    - duplicated IO is now avoided
    - added progressbars for single and multiprocessing modes

- updated __version__
parent 28235f4a
......@@ -9,7 +9,7 @@ from .components import utilities
from .components import geometry
__author__ = 'Daniel Scheffler'
__version__= '2016-11-01_01'
__version__= '2016-11-02_01'
__all__=['COREG',
'COREG_LOCAL',
......
......@@ -51,7 +51,7 @@ class imParamObj(object):
self.q = CoReg_params['q'] if not self.v else False
# set GeoArray
get_geoArr = lambda p: GeoArray(p) if not isinstance(p,GeoArray) else p
get_geoArr = lambda p: GeoArray(p) if not isinstance(p,GeoArray) else p
self.GeoArray = get_geoArr(CoReg_params['im_ref']) if imID == 'ref' else get_geoArr(CoReg_params['im_tgt'])
# set title to be used in plots
......@@ -101,7 +101,7 @@ class imParamObj(object):
#self.poly = get_footprint_polygon(self.corner_coord, fix_invalid=True) # this is the old algorithm
self.poly = self.GeoArray.footprint_poly
for XY in self.corner_coord:
assert self.poly.contains(Point(XY)) or self.poly.touches(Point(XY)), \
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('Corner coordinates of %s:\n\t%s' % (self.imName, self.corner_coord))
......@@ -112,7 +112,7 @@ class COREG(object):
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,
nodata=(None,None), calc_corners=True, multiproc=True, binary_ws=True, force_quadratic_win=True,
v=False, path_verbose_out=None, q=False, ignore_errors=False):
progress=True, v=False, path_verbose_out=None, q=False, ignore_errors=False):
"""
:param im_ref(str, GeoArray): source path (any GDAL compatible image format is supported) or GeoArray instance
of reference image
......@@ -152,6 +152,7 @@ class COREG(object):
:param multiproc(bool): enable multiprocessing (default: 1)
: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 path_verbose_out(str): an optional output directory for intermediate results
(if not given, no intermediate results are written to disk)
......@@ -195,7 +196,9 @@ class COREG(object):
self.force_quadratic_win = force_quadratic_win
self.v = v
self.path_verbose_out = path_verbose_out
self.q = q if not v else False
self.q = q if not v else False # overridden by v
self.progress = progress if not q else False # overridden by q
self.ignErr = ignore_errors
self.max_win_sz_changes = 3 # TODO: änderung der window size, falls nach max_iter kein valider match gefunden
self.ref = None # set by self.get_image_params
......@@ -245,7 +248,8 @@ class COREG(object):
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 = None if self.matchWin.boxMapYX else False
self._coreg_info = None # private attribute to be filled by self.coreg_info property
def _set_outpathes(self, im_ref, im_tgt):
......@@ -334,9 +338,9 @@ class COREG(object):
raise ImportError("This method requires the libraries 'folium' and 'geojson'. They can be installed with "
"the shell command 'pip install folium geojson'.")
refPoly = reproject_shapelyGeometry(self.ref .poly , self.ref .GeoArray.epsg, 4326)
shiftPoly = reproject_shapelyGeometry(self.shift.poly , self.shift.GeoArray.epsg, 4326)
overlapPoly = reproject_shapelyGeometry(self.overlap_poly, self.shift.GeoArray.epsg, 4326)
refPoly = reproject_shapelyGeometry(self.ref .poly , self.ref .GeoArray.epsg, 4326)
shiftPoly = reproject_shapelyGeometry(self.shift.poly , self.shift.GeoArray.epsg, 4326)
overlapPoly = reproject_shapelyGeometry(self.overlap_poly , self.shift.GeoArray.epsg, 4326)
matchWinPoly = reproject_shapelyGeometry(self.matchWin.mapPoly, self.shift.GeoArray.epsg, 4326)
m = folium.Map(location=tuple(np.array(overlapPoly.centroid.coords.xy).flatten())[::-1])
......@@ -363,8 +367,8 @@ 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 is outside of the overlap area of '\
'the two input images. Check the coordinates.'
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:
......@@ -835,17 +839,22 @@ class COREG(object):
@property
def coreg_info(self):
return {
'corrected_shifts_px' : {'x':self.x_shift_px, 'y':self.y_shift_px },
'corrected_shifts_map' : {'x':self.x_shift_map, 'y':self.y_shift_map},
'original map info' : geotransform2mapinfo(self.shift.gt, self.shift.prj),
'updated map info' : self.updated_map_info,
'reference projection' : self.ref.prj,
'reference geotransform': self.ref.gt,
'reference grid' : [ [self.ref.gt[0], self.ref.gt[0]+self.ref.gt[1]],
[self.ref.gt[3], self.ref.gt[3]+self.ref.gt[5]] ],
'reference extent' : {'cols':self.ref.xgsd, 'rows':self.ref.ygsd}, # FIXME not needed anymore
'success' : self.success}
if self._coreg_info:
return self._coreg_info
else:
self.calculate_spatial_shifts()
self._coreg_info = {
'corrected_shifts_px' : {'x':self.x_shift_px, 'y':self.y_shift_px },
'corrected_shifts_map' : {'x':self.x_shift_map, 'y':self.y_shift_map},
'original map info' : geotransform2mapinfo(self.shift.gt, self.shift.prj),
'updated map info' : self.updated_map_info,
'reference projection' : self.ref.prj,
'reference geotransform': self.ref.gt,
'reference grid' : [ [self.ref.gt[0], self.ref.gt[0]+self.ref.gt[1]],
[self.ref.gt[3], self.ref.gt[3]+self.ref.gt[5]] ],
'reference extent' : {'cols':self.ref.xgsd, 'rows':self.ref.ygsd}, # FIXME not needed anymore
'success' : self.success}
return self.coreg_info
def correct_shifts(self):
......@@ -858,6 +867,7 @@ class COREG(object):
match_gsd = self.match_gsd,
nodata = self.shift.nodata,
CPUs = None if self.mp else 1,
progress = self.progress,
v = self.v,
q = self.q)
self.deshift_results = DS.correct_shifts()
......
......@@ -25,7 +25,7 @@ class COREG_LOCAL(object):
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, v=False, q=False):
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
......@@ -54,11 +54,12 @@ class COREG_LOCAL(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: 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: 1)
:param binary_ws(bool): use binary X/Y dimensions for the matching window (default: True)
:param CPUs(int): number of CPUs to use during calculation of geometric quality grid
(default: None, which means 'all CPUs available')
:param v(bool): verbose mode (default: 0)
:param q(bool): quiet mode (default: 0)
:param progress(bool): show progress bars (default: True)
:param v(bool): verbose mode (default: False)
:param q(bool): quiet mode (default: False)
"""
self.params = dict([x for x in locals().items() if x[0] != "self" and not x[0].startswith('__')])
......@@ -81,7 +82,8 @@ class COREG_LOCAL(object):
self.bin_ws = binary_ws
self.CPUs = CPUs
self.v = v
self.q = q
self.q = q if not v else False # overridden by v
self.progress = progress if not q else False # overridden by v
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
......@@ -108,6 +110,7 @@ class COREG_LOCAL(object):
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
self.deshift_results = None # set by self.correct_shifts()
......@@ -133,7 +136,8 @@ class COREG_LOCAL(object):
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, v=self.v, q=self.q)
dir_out=self.projectDir, CPUs=self.CPUs, progress=self.progress,
v=self.v, q=self.q)
return self._quality_grid
......@@ -218,9 +222,11 @@ class COREG_LOCAL(object):
center_lon, center_lat = (lon_min+lon_max)/2, (lat_min+lat_max)/2
# get image to plot
image2plot = self.im2shift[:]
image2plot = self.im2shift[0] # FIXME hardcoded band
from py_tools_ds.ptds.geo.raster.reproject import warp_ndarray
image2plot, gt, prj = warp_ndarray(image2plot, self.im2shift.geotransform, self.im2shift.projection, in_nodata=self.nodata[1], out_nodata=self.nodata[1],
image2plot, gt, prj = warp_ndarray(image2plot, self.im2shift.geotransform, self.im2shift.projection,
in_nodata=self.nodata[1], out_nodata=self.nodata[1],
out_XYdims=(1000, 1000), q=True, out_prj='epsg:3857') # image must be transformed into web mercator projection
# create map
......@@ -243,15 +249,32 @@ class COREG_LOCAL(object):
return map_osm
def correct_shifts(self, max_GCP_count=None):
@property
def coreg_info(self):
if self._coreg_info:
return self._coreg_info
else:
self._coreg_info = {
'GCPList' : self.quality_grid.GCPList,
'reference projection' : self.imref.prj,
'reference geotransform': self.im2shift.gt,
'reference grid' : [ [self.imref.gt[0], self.imref.gt[0]+self.imref.gt[1]],
[self.imref.gt[3], self.imref.gt[3]+self.imref.gt[5]] ],
'reference extent' : {'cols':self.imref.xgsd, 'rows':self.imref.ygsd}, # FIXME not needed anymore
#'success' : self.success
}
return self.coreg_info
def correct_shifts(self, max_GCP_count=None, cliptoextent=False):
"""Performs a local shift correction using all points from the previously calculated geometric quality grid
that contain valid matches as GCP points.
:param max_GCP_count: <int> maximum number of GCPs to use
:param cliptoextent: <bool> whether to clip the output image to its real extent
:return:
"""
coreg_info = self.COREG_obj.coreg_info
coreg_info['GCPList'] = self.quality_grid.GCPList
coreg_info = self.coreg_info
if max_GCP_count:
coreg_info['GCPList'] = coreg_info['GCPList'][:max_GCP_count] # TODO should be a random sample
......@@ -261,10 +284,9 @@ class COREG_LOCAL(object):
fmt_out = self.fmt_out,
out_gsd = (self.im2shift.xgsd,self.im2shift.ygsd),
align_grids = True,
#cliptoextent = True, # why?
#cliptoextent = False,
cliptoextent = cliptoextent,
#clipextent = self.im2shift.box.boxMapYX,
#options = '-wm 10000 -order 1',
progress = self.progress,
v = self.v,
q = self.q)
......
......@@ -2,8 +2,6 @@
__author__='Daniel Scheffler'
import collections
import os
import tempfile
import time
import warnings
......@@ -12,8 +10,6 @@ try:
import gdal
except ImportError:
from osgeo import gdal
import numpy as np
import rasterio
# internal modules
from py_tools_ds.ptds import GeoArray
......@@ -22,7 +18,6 @@ from py_tools_ds.ptds.geo.coord_grid import is_coord_grid_equal
from py_tools_ds.ptds.geo.projection import prj_equal
from py_tools_ds.ptds.geo.raster.reproject import warp_ndarray
from py_tools_ds.ptds.numeric.vector import find_nearest
from py_tools_ds.ptds.processing.shell import subcall_with_output
_dict_rspAlg_rsp_Int = {'nearest': 0, 'bilinear': 1, 'cubic': 2, 'cubic_spline': 3, 'lanczos': 4, 'average': 5,
'mode': 6, 'max': 7, 'min': 8 , 'med': 9, 'q1':10, 'q2':11}
......@@ -33,7 +28,8 @@ class DESHIFTER(object):
Deshift an image array or one of its products by applying the coregistration info calculated by COREG class.
:param im2shift: <path,GeoArray> path of an image to be de-shifted or alternatively a GeoArray object
:param coreg_results: <dict> the results of the co-registration as given by COREG.coreg_info
:param coreg_results: <dict> the results of the co-registration as given by COREG.coreg_info or
COREG_LOCAL.coreg_info respectively
:Keyword Arguments:
- path_out(str): /output/directory/filename for coregistered results
......@@ -53,13 +49,12 @@ class DESHIFTER(object):
- 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)
- warp_alg(str): 'GDAL_cmd' or 'GDAL_lib' (default = 'GDAL_lib')
- cliptoextent (bool): True: clip the input image to its actual bounds while deleting possible no data
areas outside of the actual bounds, default = True
- clipextent (list): xmin, ymin, xmax, ymax - if given the calculation of the actual bounds is skipped.
The given coordinates are automatically snapped to the output grid.
- tempDir(str): directory to be used for tempfiles (default: /dev/shm/)
- CPUs(int): number of CPUs to use (default: None, which means 'all CPUs available')
- progress(bool): show progress bars (default: True)
- v(bool): verbose mode (default: False)
- q(bool): quiet mode (default: False)
......@@ -69,10 +64,11 @@ class DESHIFTER(object):
self.shift_prj = self.im2shift.projection
self.shift_gt = list(self.im2shift.geotransform)
self.GCPList = coreg_results['GCPList'] if 'GCPList' in coreg_results else None
mapI = coreg_results['updated map info']
self.updated_map_info = mapI if mapI else geotransform2mapinfo(self.shift_gt, self.shift_prj)
self.original_map_info = coreg_results['original map info']
self.updated_gt = mapinfo2geotransform(self.updated_map_info) if mapI else self.shift_gt
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)
self.original_map_info = coreg_results['original map info']
self.updated_gt = mapinfo2geotransform(self.updated_map_info) if mapI else self.shift_gt
self.ref_gt = coreg_results['reference geotransform']
self.ref_grid = coreg_results['reference grid']
self.ref_prj = coreg_results['reference projection']
......@@ -84,23 +80,19 @@ class DESHIFTER(object):
self.band2process = kwargs.get('band2process', None) # starts with 1 # FIXME warum?
self.nodata = kwargs.get('nodata' , self.im2shift.nodata)
self.align_grids = kwargs.get('align_grids' , False)
tempAsENVI = kwargs.get('tempAsENVI' , False)
self.outFmt = 'VRT' if not tempAsENVI else 'ENVI' # FIXME eliminate that
self.rspAlg = kwargs.get('resamp_alg' , 'cubic')
self.warpAlg = kwargs.get('warp_alg' , 'GDAL_lib')
self.cliptoextent = kwargs.get('cliptoextent', True)
self.clipextent = kwargs.get('clipextent' , None)
self.tempDir = kwargs.get('tempDir' , '/dev/shm/')
self.CPUs = kwargs.get('CPUs' , None)
self.v = kwargs.get('v' , False)
self.q = kwargs.get('q' , False) if not self.v else False
self.q = kwargs.get('q' , False) if not self.v else False # overridden by v
self.progress = kwargs.get('progress' , True) if not self.q else False # overridden by q
self.out_grid = self._get_out_grid(kwargs) # 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
assert self.warpAlg in ['GDAL_cmd', 'GDAL_lib']
# set defaults for general class attributes
self.is_shifted = False # this is not included in COREG.coreg_info
......@@ -193,6 +185,9 @@ class DESHIFTER(object):
def correct_shifts(self):
# type: (DESHIFTER) -> collections.OrderedDict
if not self.q:
print('Correcting geometric shifts...')
t_start = time.time()
equal_prj = prj_equal(self.ref_prj,self.shift_prj)
......@@ -225,77 +220,83 @@ class DESHIFTER(object):
else: # FIXME equal_prj==False ist noch NICHT implementiert
"""RESAMPLING NEEDED"""
if self.warpAlg=='GDAL_cmd':
warnings.warn('This method has not been tested in its current state!')
# FIXME nicht multiprocessing-fähig, weil immer kompletter array gewarpt wird und sich ergebnisse gegenseitig überschreiben
# create tempfile
fd, path_tmp = tempfile.mkstemp(prefix='CoReg_Sat', suffix=self.outFmt, dir=self.tempDir)
os.close(fd)
t_extent = " -te %s %s %s %s" %self._get_out_extent()
xgsd, ygsd = self.out_gsd
cmd = "gdalwarp -r %s -tr %s %s -t_srs '%s' -of %s %s %s -srcnodata %s -dstnodata %s -overwrite%s"\
%(self.rspAlg, xgsd,ygsd,self.ref_prj,self.outFmt,self.im2shift.filePath,
path_tmp, self.nodata, self.nodata, t_extent)
out, exitcode, err = subcall_with_output(cmd)
if exitcode!=1 and os.path.exists(path_tmp):
"""update map info, arr_shifted, geotransform and projection"""
ds_shifted = gdal.OpenShared(path_tmp) if self.outFmt == 'VRT' else gdal.Open(path_tmp)
self.shift_gt, self.shift_prj = ds_shifted.GetGeoTransform(), ds_shifted.GetProjection()
self.updated_map_info = geotransform2mapinfo(self.shift_gt,self.shift_prj)
print('reading from', ds_shifted.GetDescription())
if self.band2process is None:
dim2RowsColsBands = lambda A: np.swapaxes(np.swapaxes(A,0,2),0,1) # rasterio.open(): [bands,rows,cols]
self.arr_shifted = dim2RowsColsBands(rasterio.open(path_tmp).read())
else:
self.arr_shifted = rasterio.open(path_tmp).read(self.band2process)
self.GeoArray_shifted = GeoArray(self.arr_shifted,tuple(self.shift_gt), self.shift_prj)
self.is_shifted = True
self.is_resampled = True
ds_shifted = None
[gdal.Unlink(p) for p in [path_tmp] if os.path.exists(p)] # delete tempfiles
else:
print("\n%s\nCommand was: '%s'" %(err.decode('utf8'),cmd))
[gdal.Unlink(p) for p in [path_tmp] if os.path.exists(p)] # delete tempfiles
self.tracked_errors.append(RuntimeError('Resampling failed.'))
raise self.tracked_errors[-1]
# TODO implement output writer
elif self.warpAlg=='GDAL_lib':
# if self.warpAlg=='GDAL_cmd':
# warnings.warn('This method has not been tested in its current state!')
# # FIX ME nicht multiprocessing-fähig, weil immer kompletter array gewarpt wird und sich ergebnisse gegenseitig überschreiben
# # create tempfile
# fd, path_tmp = tempfile.mkstemp(prefix='CoReg_Sat', suffix=self.outFmt, dir=self.tempDir)
# os.close(fd)
#
# t_extent = " -te %s %s %s %s" %self._get_out_extent()
# xgsd, ygsd = self.out_gsd
# cmd = "gdalwarp -r %s -tr %s %s -t_srs '%s' -of %s %s %s -srcnodata %s -dstnodata %s -overwrite%s"\
# %(self.rspAlg, xgsd,ygsd,self.ref_prj,self.outFmt,self.im2shift.filePath,
# path_tmp, self.nodata, self.nodata, t_extent)
# out, exitcode, err = subcall_with_output(cmd)
#
# if exitcode!=1 and os.path.exists(path_tmp):
# """update map info, arr_shifted, geotransform and projection"""
# ds_shifted = gdal.OpenShared(path_tmp) if self.outFmt == 'VRT' else gdal.Open(path_tmp)
# self.shift_gt, self.shift_prj = ds_shifted.GetGeoTransform(), ds_shifted.GetProjection()
# self.updated_map_info = geotransform2mapinfo(self.shift_gt,self.shift_prj)
#
# print('reading from', ds_shifted.GetDescription())
# if self.band2process is None:
# dim2RowsColsBands = lambda A: np.swapaxes(np.swapaxes(A,0,2),0,1) # rasterio.open(): [bands,rows,cols]
# self.arr_shifted = dim2RowsColsBands(rasterio.open(path_tmp).read())
# else:
# self.arr_shifted = rasterio.open(path_tmp).read(self.band2process)
#
# self.GeoArray_shifted = GeoArray(self.arr_shifted,tuple(self.shift_gt), self.shift_prj)
# self.is_shifted = True
# self.is_resampled = True
#
# ds_shifted = None
# [gdal.Unlink(p) for p in [path_tmp] if os.path.exists(p)] # delete tempfiles
# else:
# print("\n%s\nCommand was: '%s'" %(err.decode('utf8'),cmd))
# [gdal.Unlink(p) for p in [path_tmp] if os.path.exists(p)] # delete tempfiles
# self.tracked_errors.append(RuntimeError('Resampling failed.'))
# raise self.tracked_errors[-1]
#
# # TO DO implement output writer
in_arr = self.im2shift[self.band2process] if self.band2process else self.im2shift[:]
if not self.GCPList:
# apply XY-shifts to shift_gt
in_arr = self.im2shift[self.band2process] if self.band2process else self.im2shift[:]
if not self.GCPList:
self.shift_gt[0], self.shift_gt[3] = self.updated_gt[0], self.updated_gt[3]
# get resampled array
out_arr, out_gt, out_prj = \
warp_ndarray(in_arr, self.shift_gt, self.shift_prj, self.ref_prj,
rspAlg = _dict_rspAlg_rsp_Int[self.rspAlg],
in_nodata = self.nodata,
out_nodata = self.nodata,
out_gsd = self.out_gsd,
out_bounds = self._get_out_extent(),
gcpList = self.GCPList,
polynomialOrder = None,
options = None, #'-refine_gcps 500',
CPUs = self.CPUs,
q = self.q)
self.updated_projection = out_prj
self.arr_shifted = out_arr
self.updated_map_info = geotransform2mapinfo(out_gt,out_prj)
self.shift_gt = mapinfo2geotransform(self.updated_map_info)
self.GeoArray_shifted = GeoArray(self.arr_shifted, tuple(self.shift_gt), self.updated_projection)
self.is_shifted = True
self.is_resampled = True
if self.path_out:
GeoArray(out_arr, out_gt, out_prj).save(self.path_out,fmt=self.fmt_out)
self.shift_gt[0], self.shift_gt[3] = self.updated_gt[0], self.updated_gt[3]
# get resampled array
out_arr, out_gt, out_prj = \
warp_ndarray(in_arr, self.shift_gt, self.shift_prj, self.ref_prj,
rspAlg = _dict_rspAlg_rsp_Int[self.rspAlg],
in_nodata = self.nodata,
out_nodata = self.nodata,
out_gsd = self.out_gsd,
out_bounds = self._get_out_extent(),
gcpList = self.GCPList,
#polynomialOrder = str(3),
#options = '-refine_gcps 500 1.9',
#warpOptions= ['-refine_gcps 500 1.9'],
#options = '-wm 10000',# -order 3',
#options = ['-order 3'],
# options = ['GDAL_CACHEMAX 800 '],
#warpMemoryLimit=125829120, # 120MB
CPUs = self.CPUs,
progress = self.progress,
q = self.q)
self.updated_projection = out_prj
self.arr_shifted = out_arr
self.updated_map_info = geotransform2mapinfo(out_gt,out_prj)
self.shift_gt = mapinfo2geotransform(self.updated_map_info)
self.GeoArray_shifted = GeoArray(self.arr_shifted, tuple(self.shift_gt), self.updated_projection)
self.is_shifted = True
self.is_resampled = True
if self.path_out:
GeoArray(out_arr, out_gt, out_prj).save(self.path_out,fmt=self.fmt_out)
if self.v: print('Time for shift correction: %.2fs' %(time.time()-t_start))
return self.deshift_results
......
......@@ -5,6 +5,7 @@ import collections
import multiprocessing
import os
import warnings
import time
# custom
try:
......@@ -18,10 +19,10 @@ from shapely.geometry import Point
# internal modules
from .CoReg import COREG
from . import io as IO
from py_tools_ds.ptds import GeoArray
from py_tools_ds.ptds.geo.projection import isProjectedOrGeographic, get_UTMzone
from py_tools_ds.ptds.io.pathgen import get_generic_outpath
from . import io as IO
from py_tools_ds.ptds.geo.projection import isProjectedOrGeographic, get_UTMzone
from py_tools_ds.ptds.io.pathgen import get_generic_outpath
from py_tools_ds.ptds.processing.progress_mon import printProgress
......@@ -30,7 +31,7 @@ global_shared_im2shift = None
class Geom_Quality_Grid(object):
def __init__(self, COREG_obj, grid_res, outFillVal=-9999, dir_out=None, CPUs=None, v=False, q=False):
def __init__(self, COREG_obj, grid_res, outFillVal=-9999, 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
......@@ -45,6 +46,7 @@ class Geom_Quality_Grid(object):
to the individual methods
: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)
:param v(bool): verbose mode (default: 0)
:param q(bool): quiet mode (default: 0)
"""
......@@ -57,7 +59,8 @@ class Geom_Quality_Grid(object):
self.outFillVal = outFillVal
self.CPUs = CPUs
self.v = v
self.q = q
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
......@@ -125,7 +128,8 @@ class Geom_Quality_Grid(object):
# imX, imY = mapXY2imXY(coreg_kwargs['wp'], im.gt)
# if im.GeoArray[int(imY), int(imX), im.band4match]==im.nodata,\
# return
assert global_shared_imref is not None
assert global_shared_im2shift is not None
CR = COREG(global_shared_imref, global_shared_im2shift, multiproc=False, **coreg_kwargs)
CR.calculate_spatial_shifts()
......@@ -183,10 +187,12 @@ 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.ref[:,:,self.COREG_obj.ref.band4match-1], self.ref.geotransform, self.ref.projection)
global_shared_im2shift = \
GeoArray(self.shift[:,:,self.COREG_obj.shift.band4match-1], self.shift.geotransform,self.shift.projection)
assert self.ref .footprint_poly # this also checks for mask_nodata and nodata value
assert self.shift.footprint_poly
if not self.ref .is_inmem: self.ref.cache_array_subset(self.ref [self.COREG_obj.ref .band4match])
if not self.shift.is_inmem: self.ref.cache_array_subset(self.shift[self.COREG_obj.shift.band4match])
global_shared_imref = self.ref
global_shared_im2shift = self.shift
# get all variations of kwargs for coregistration
get_coreg_kwargs = lambda pID, wp: {
......@@ -211,21 +217,34 @@ class Geom_Quality_Grid(object):
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))
t0 = time.time()
with multiprocessing.Pool(self.CPUs) as pool:
results = pool.map(self._get_spatial_shifts, list_coreg_kwargs)
if self.q or not self.progress:
results = pool.map(self._get_spatial_shifts, list_coreg_kwargs)
else:
results = pool.map_async(self._get_spatial_shifts, list_coreg_kwargs, chunksize=1)
while True:
time.sleep(.1)
numberDone = len(GDF)-results._number_left # this does not really represent the remaining tasks but the remaining chunks -> thus chunksize=1
printProgress(percent=numberDone/len(GDF)*100, barLength=50, prefix='\tprogress:',
suffix='[%s/%s] Complete %.2f sek'%(numberDone,len(GDF), time.time()-t0))
if results.ready():
results = results.get() # FIXME in some cases the code hangs here ==> x
break
else:
if not self.q:
print("Calculating geometric quality grid (%s points) in mode 'singleprocessing'..." %len(GDF))
results = np.empty((len(geomPoints),9))
for i,coreg_kwargs in enumerate(list_coreg_kwargs):
#print(argset[1])