Commit 02e7bb36 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Revised error handling and added additional check for projection. Some...

Revised error handling and added additional check for projection. Some PEP8-editing. Test artifacts are now always saved.
parent f64a77fe
Pipeline #1023 failed with stages
in 6 minutes and 52 seconds
......@@ -32,6 +32,7 @@ test_arosics:
- docs/_build/html/
- nosetests.html
- nosetests.xml
when: always
test_arosics_install:
......
# -*- coding: utf-8 -*-
__author__='Daniel Scheffler'
import os
import re
......@@ -15,6 +14,7 @@ try:
except ImportError:
from osgeo import gdal
import numpy as np
try:
import pyfftw
except ImportError:
......@@ -24,9 +24,9 @@ from skimage.exposure import rescale_intensity
# internal modules
from .DeShifter import DESHIFTER, _dict_rspAlg_rsp_Int
from . import geometry as GEO
from . import io as IO
from . import plotting as PLT
from . import geometry as GEO
from . import io as IO
from . import plotting as PLT
from geoarray import GeoArray
from py_tools_ds.convenience.object_oriented import alias_property
......@@ -41,7 +41,7 @@ from py_tools_ds.geo.map_info import geotransform2mapinfo
from py_tools_ds.numeric.vector import find_nearest
from py_tools_ds.similarity.raster import calc_ssim
__author__ = 'Daniel Scheffler'
class GeoArray_CoReg(GeoArray):
......@@ -52,20 +52,20 @@ class GeoArray_CoReg(GeoArray):
# run GeoArray init
path_or_geoArr = CoReg_params['im_ref'] if imID == 'ref' else CoReg_params['im_tgt']
nodata = CoReg_params['nodata'][0 if imID == 'ref' else 1]
progress = CoReg_params['progress']
q = CoReg_params['q'] if not CoReg_params['v'] else False
nodata = CoReg_params['nodata'][0 if imID == 'ref' else 1]
progress = CoReg_params['progress']
q = CoReg_params['q'] if not CoReg_params['v'] else False
super(GeoArray_CoReg, self).__init__(path_or_geoArr, nodata=nodata, progress=progress, q=q)
self.imID = imID
self.imID = imID
self.imName = 'reference image' if imID == 'ref' else 'image to be shifted'
self.v = CoReg_params['v']
self.v = CoReg_params['v']
assert isinstance(self, GeoArray), \
'Something went wrong with the creation of GeoArray instance for the %s. The created ' \
'instance does not seem to belong to the GeoArray class. If you are working in Jupyter Notebook, reset ' \
'the kernel and try again.' %self.imName
'the kernel and try again.' % self.imName
# set title to be used in plots
self.title = os.path.basename(self.filePath) if self.filePath else self.imName
......@@ -76,14 +76,15 @@ class GeoArray_CoReg(GeoArray):
assert self.gt, 'The %s has no map information.' % self.imName
# set band4match
self.band4match = (CoReg_params['r_b4match'] if imID == 'ref' else CoReg_params['s_b4match'])-1
assert self.bands >= self.band4match+1 >= 1, "The %s has %s %s. So its band number to match must be %s%s. " \
"Got %s." % (self.imName, self.bands, 'bands' if self.bands > 1 else 'band', 'between 1 and '
if self.bands > 1 else '', self.bands, self.band4match)
self.band4match = (CoReg_params['r_b4match'] if imID == 'ref' else CoReg_params['s_b4match']) - 1
assert self.bands >= self.band4match + 1 >= 1, \
"The %s has %s %s. So its band number to match must be %s%s. Got %s." \
% (self.imName, self.bands, 'bands' if self.bands > 1 else
'band', 'between 1 and ' if self.bands > 1 else '', self.bands, self.band4match)
# 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')]
given_corner_coord = CoReg_params['data_corners_%s' % ('ref' if imID == 'ref' else 'tgt')]
if given_footprint_poly:
self.footprint_poly = given_footprint_poly
......@@ -91,7 +92,7 @@ class GeoArray_CoReg(GeoArray):
self.footprint_poly = Polygon(given_corner_coord)
elif not CoReg_params['calc_corners']:
# use the image extent
self.footprint_poly = Polygon(get_corner_coordinates(gt=self.gt, cols=self.cols,rows=self.rows))
self.footprint_poly = Polygon(get_corner_coordinates(gt=self.gt, cols=self.cols, rows=self.rows))
else:
# footprint_poly is calculated automatically by GeoArray
if not CoReg_params['q']:
......@@ -104,21 +105,19 @@ class GeoArray_CoReg(GeoArray):
# add bad data mask
given_mask = CoReg_params['mask_baddata_%s' % ('ref' if imID == 'ref' else 'tgt')]
if given_mask:
self.mask_baddata = given_mask # runs GeoArray.mask_baddata.setter -> sets it to BadDataMask()
poly = alias_property('footprint_poly') # ensures that self.poly is updated if self.footprint_poly is updated
self.mask_baddata = given_mask # runs GeoArray.mask_baddata.setter -> sets it to BadDataMask()
poly = alias_property('footprint_poly') # ensures that self.poly is updated if self.footprint_poly is updated
class COREG(object):
"""See help(COREG) for documentation!"""
def __init__(self, im_ref, im_tgt, path_out=None, fmt_out='ENVI', out_crea_options=None, r_b4match=1, s_b4match=1,
wp=(None,None), ws=(256, 256), max_iter=5, max_shift=5, align_grids=False, match_gsd=False,
wp=(None, None), ws=(256, 256), 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, binary_ws=True, mask_baddata_ref=None, mask_baddata_tgt=None,
nodata=(None, None), calc_corners=True, binary_ws=True, mask_baddata_ref=None, mask_baddata_tgt=None,
CPUs=None, force_quadratic_win=True, progress=True, v=False, path_verbose_out=None, q=False,
ignore_errors=False):
......@@ -144,25 +143,28 @@ class COREG(object):
:param ws(tuple): custom matching window size [pixels] (default: (256,256))
: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 align_grids(bool): align the coordinate grids of the image to be and the reference image (default: 0)
:param align_grids(bool): align the coordinate grids of the image to be and the reference image
(default: 0)
:param match_gsd(bool): match the output pixel size to pixel size of the reference image (default: 0)
:param out_gsd(tuple): xgsd ygsd: set the output pixel size in map units
(default: original pixel size of the image to be shifted)
: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
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)
(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),
: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)
: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.
......@@ -199,18 +201,22 @@ class COREG(object):
is None
"""
self.params = dict([x for x in locals().items() if x[0] != "self"])
self.params = dict([x for x in locals().items() if x[0] != "self"])
# assertions
assert gdal.GetDriverByName(fmt_out), "'%s' is not a supported GDAL driver." % 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.'
if data_corners_ref and not isinstance(data_corners_ref[0], list): # group if not [[x,y],[x,y]..] but [x,y,x,y,]
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_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]..]
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))
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]:
assert rspAlg in _dict_rspAlg_rsp_Int.keys(), "'%s' is not a supported resampling algorithm." % rspAlg
if resamp_alg_calc in ['average', 5] and (v or not q):
......@@ -218,88 +224,90 @@ class COREG(object):
"affect the precision of the calculated spatial shifts! It is highly recommended to "
"choose another resampling algorithm.")
self.path_out = path_out # updated by self.set_outpathes
self.fmt_out = fmt_out
self.out_creaOpt = out_crea_options
self.win_pos_XY = wp # updated by self.get_opt_winpos_winsize()
self.win_size_XY = ws # updated by self.get_opt_winpos_winsize()
self.max_iter = max_iter
self.max_shift = max_shift
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.path_out = path_out # updated by self.set_outpathes
self.fmt_out = fmt_out
self.out_creaOpt = out_crea_options
self.win_pos_XY = wp # updated by self.get_opt_winpos_winsize()
self.win_size_XY = ws # updated by self.get_opt_winpos_winsize()
self.max_iter = max_iter
self.max_shift = max_shift
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 \
if isinstance(resamp_alg_deshift, str) else _dict_rspAlg_rsp_Int[resamp_alg_deshift]
self.rspAlg_calc = resamp_alg_calc \
self.rspAlg_calc = resamp_alg_calc \
if isinstance(resamp_alg_calc, str) else _dict_rspAlg_rsp_Int[resamp_alg_calc]
self.calc_corners = calc_corners
self.CPUs = CPUs
self.bin_ws = binary_ws
self.calc_corners = calc_corners
self.CPUs = CPUs
self.bin_ws = binary_ws
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 # 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
self.shift = None # set by self.get_image_params
self.matchBox = None # set by self.get_clip_window_properties() => boxObj
self.otherBox = None # set by self.get_clip_window_properties() => boxObj
self.matchWin = None # set by self._get_image_windows_to_match() => GeoArray
self.otherWin = None # set by self._get_image_windows_to_match() => GeoArray
self.overlap_poly = None # set by self._get_overlap_properties()
self.overlap_percentage = None # set by self._get_overlap_properties()
self.overlap_area = None # set by self._get_overlap_properties()
self.imfft_gsd = None # set by self.get_clip_window_properties()
self.fftw_works = None # set by self._calc_shifted_cross_power_spectrum()
self.fftw_win_size_YX = None # set by calc_shifted_cross_power_spectrum()
self.x_shift_px = None # always in shift image units (image coords) # set by calculate_spatial_shifts()
self.y_shift_px = None # always in shift image units (image coords) # set by calculate_spatial_shifts()
self.x_shift_map = None # set by self.get_updated_map_info()
self.y_shift_map = None # set by self.get_updated_map_info()
self.vec_length_map = None
self.vec_angle_deg = None
self.updated_map_info = None # set by self.get_updated_map_info()
self.ssim_orig = None # set by self._validate_ssim_improvement()
self.ssim_deshifted = None # set by self._validate_ssim_improvement()
self._ssim_improved = None # private attribute to be filled by self.ssim_improved
self.shift_reliability = None # set by self.calculate_spatial_shifts()
self.tracked_errors = [] # expanded each time an error occurs
self.success = None # default
self.deshift_results = None # set by self.correct_shifts()
self.v = v
self.path_verbose_out = path_verbose_out
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
self.shift = None # set by self.get_image_params
self.matchBox = None # set by self.get_clip_window_properties() => boxObj
self.otherBox = None # set by self.get_clip_window_properties() => boxObj
self.matchWin = None # set by self._get_image_windows_to_match() => GeoArray
self.otherWin = None # set by self._get_image_windows_to_match() => GeoArray
self.overlap_poly = None # set by self._get_overlap_properties()
self.overlap_percentage = None # set by self._get_overlap_properties()
self.overlap_area = None # set by self._get_overlap_properties()
self.imfft_gsd = None # set by self.get_clip_window_properties()
self.fftw_works = None # set by self._calc_shifted_cross_power_spectrum()
self.fftw_win_size_YX = None # set by calc_shifted_cross_power_spectrum()
self.x_shift_px = None # always in shift image units (image coords) # set by calculate_spatial_shifts()
self.y_shift_px = None # always in shift image units (image coords) # set by calculate_spatial_shifts()
self.x_shift_map = None # set by self.get_updated_map_info()
self.y_shift_map = None # set by self.get_updated_map_info()
self.vec_length_map = None
self.vec_angle_deg = None
self.updated_map_info = None # set by self.get_updated_map_info()
self.ssim_orig = None # set by self._validate_ssim_improvement()
self.ssim_deshifted = None # set by self._validate_ssim_improvement()
self._ssim_improved = None # private attribute to be filled by self.ssim_improved
self.shift_reliability = None # set by self.calculate_spatial_shifts()
self.tracked_errors = [] # expanded each time an error occurs
self.success = None # default
self.deshift_results = None # set by self.correct_shifts()
gdal.AllRegister()
self._get_image_params()
self._set_outpathes(im_ref, im_tgt)
self.grid2use = 'ref' if self.shift.xgsd <= self.ref.xgsd else 'shift'
if self.v: print('resolutions: ', self.ref.xgsd, self.shift.xgsd)
if self.v:
print('resolutions: ', self.ref.xgsd, self.shift.xgsd)
self._get_overlap_properties()
if self.v and self.path_verbose_out:
IO.write_shp(os.path.join(self.path_verbose_out, 'poly_imref.shp'), self.ref.poly, self.ref.prj)
IO.write_shp(os.path.join(self.path_verbose_out, 'poly_im2shift.shp'), self.shift.poly, self.shift.prj)
IO.write_shp(os.path.join(self.path_verbose_out, 'overlap_poly.shp'), self.overlap_poly, self.ref.prj)
IO.write_shp(os.path.join(self.path_verbose_out, 'poly_imref.shp'), self.ref.poly, self.ref.prj)
IO.write_shp(os.path.join(self.path_verbose_out, 'poly_im2shift.shp'), self.shift.poly, self.shift.prj)
IO.write_shp(os.path.join(self.path_verbose_out, 'overlap_poly.shp'), self.overlap_poly, self.ref.prj)
### FIXME: transform_mapPt1_to_mapPt2(im2shift_center_map, ds_imref.GetProjection(), ds_im2shift.GetProjection()) # später basteln für den fall, dass projektionen nicht gleich sind
# FIXME: transform_mapPt1_to_mapPt2(im2shift_center_map, ds_imref.GetProjection(), ds_im2shift.GetProjection())
# FIXME später basteln für den fall, dass projektionen nicht gleich sind
# get_clip_window_properties
self._get_opt_winpos_winsize()
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() # sets self.matchBox, self.otherBox and much more
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() # sets self.matchBox, self.otherBox and much more
if self.v and self.path_verbose_out and self.matchBox.mapPoly and self.success is not False:
IO.write_shp(os.path.join(self.path_verbose_out, 'poly_matchWin.shp'),
self.matchBox.mapPoly, self.matchBox.prj)
self.success = False if self.success is False or not self.matchBox.boxMapYX else None
self._coreg_info = None # private attribute to be filled by self.coreg_info property
self.success = False if self.success is False or not self.matchBox.boxMapYX else None
self._coreg_info = None # private attribute to be filled by self.coreg_info property
def _handle_error(self, error, warn=False, warnMsg=None):
"""Appends the given error to self.tracked_errors, sets self.success to False and raises the error in case
......@@ -318,16 +326,15 @@ class COREG(object):
if self.ignErr and warn:
warnMsg = repr(error) if not warnMsg else warnMsg
print('\nWARNING: '+warnMsg)
print('\nWARNING: ' + warnMsg)
if not self.ignErr:
raise error
def _set_outpathes(self, im_ref, im_tgt):
assert isinstance(im_ref, (GeoArray, str)) and isinstance(im_tgt, (GeoArray, str)),\
assert isinstance(im_ref, (GeoArray, str)) and isinstance(im_tgt, (GeoArray, str)), \
'COREG._set_outpathes() expects two file pathes (string) or two instances of the ' \
'GeoArray class. Received %s and %s.' %(type(im_ref), type(im_tgt))
'GeoArray class. Received %s and %s.' % (type(im_ref), type(im_tgt))
get_baseN = lambda path: os.path.splitext(os.path.basename(path))[0]
......@@ -335,7 +342,7 @@ class COREG(object):
path_im_ref = im_ref.filePath if isinstance(im_ref, GeoArray) else im_ref
path_im_tgt = im_tgt.filePath if isinstance(im_tgt, GeoArray) else im_tgt
if self.path_out: # this also applies to self.path_out='auto'
if self.path_out: # this also applies to self.path_out='auto'
if self.path_out == 'auto':
dir_out, fName_out = os.path.dirname(path_im_tgt), ''
......@@ -355,13 +362,13 @@ class COREG(object):
dir_out = os.path.dirname(path_im_ref)
if not fName_out:
ext = 'bsq' if self.fmt_out=='ENVI' else \
gdal.GetDriverByName(self.fmt_out).GetMetadataItem(gdal.DMD_EXTENSION)
fName_out = fName_out if not fName_out in ['.',''] else '%s__shifted_to__%s' \
%(get_baseN(path_im_tgt), get_baseN(path_im_ref))
fName_out = fName_out+'.%s'%ext if ext else fName_out
ext = 'bsq' if self.fmt_out == 'ENVI' else \
gdal.GetDriverByName(self.fmt_out).GetMetadataItem(gdal.DMD_EXTENSION)
fName_out = fName_out if fName_out not in ['.', ''] else \
'%s__shifted_to__%s' % (get_baseN(path_im_tgt), get_baseN(path_im_ref))
fName_out = fName_out + '.%s' % ext if ext else fName_out
self.path_out = os.path.abspath(os.path.join(dir_out,fName_out))
self.path_out = os.path.abspath(os.path.join(dir_out, fName_out))
assert ' ' not in self.path_out, \
"The path of the output image contains whitespaces. This is not supported by GDAL."
......@@ -378,8 +385,9 @@ class COREG(object):
if self.path_out:
self.path_verbose_out = os.path.dirname(self.path_out)
else:
self.path_verbose_out = os.path.abspath(os.path.join(os.path.curdir,
'CoReg_verboseOut__%s__shifted_to__%s' % (get_baseN(path_im_tgt), get_baseN(path_im_ref))))
self.path_verbose_out = \
os.path.abspath(os.path.join(os.path.curdir, 'CoReg_verboseOut__%s__shifted_to__%s'
% (get_baseN(path_im_tgt), get_baseN(path_im_ref))))
elif dirname_out and not dir_out:
self.path_verbose_out = os.path.abspath(os.path.join(os.path.curdir, dirname_out))
......@@ -389,51 +397,52 @@ class COREG(object):
else:
self.path_verbose_out = None
if self.path_verbose_out and not os.path.isdir(self.path_verbose_out): os.makedirs(self.path_verbose_out)
if self.path_verbose_out and not os.path.isdir(self.path_verbose_out):
os.makedirs(self.path_verbose_out)
def _get_image_params(self):
self.ref = GeoArray_CoReg(self.params,'ref')
self.shift = GeoArray_CoReg(self.params,'shift')
self.ref = GeoArray_CoReg(self.params, 'ref')
self.shift = GeoArray_CoReg(self.params, 'shift')
assert self.ref.prj, 'The reference image has no projection.'
assert self.shift.prj, 'The target image has no projection.'
assert prj_equal(self.ref.prj, self.shift.prj), \
'Input projections are not equal. Different projections are currently not supported. Got %s / %s.'\
%(get_proj4info(proj=self.ref.prj), get_proj4info(proj=self.shift.prj))
'Input projections are not equal. Different projections are currently not supported. Got %s / %s.' \
% (get_proj4info(proj=self.ref.prj), get_proj4info(proj=self.shift.prj))
def _get_overlap_properties(self):
overlap_tmp = get_overlap_polygon(self.ref.poly, self.shift.poly, self.v)
self.overlap_poly = overlap_tmp['overlap poly'] # has to be in reference projection
self.overlap_percentage = overlap_tmp['overlap percentage']
self.overlap_area = overlap_tmp['overlap area']
overlap_tmp = get_overlap_polygon(self.ref.poly, self.shift.poly, self.v)
self.overlap_poly = overlap_tmp['overlap poly'] # has to be in reference projection
self.overlap_percentage = overlap_tmp['overlap percentage']
self.overlap_area = overlap_tmp['overlap area']
assert self.overlap_poly, 'The input images have no spatial overlap.'
# overlap are must at least cover 16*16 pixels
px_area = self.ref.xgsd * self.ref.ygsd if self.grid2use=='ref' else self.shift.xgsd * self.shift.ygsd
px_covered = self.overlap_area/px_area
assert px_covered > 16*16, 'Overlap area covers only %s pixels. At least 16*16 pixels are needed.' %px_covered
px_area = self.ref.xgsd * self.ref.ygsd if self.grid2use == 'ref' else self.shift.xgsd * self.shift.ygsd
px_covered = self.overlap_area / px_area
assert px_covered > 16 * 16, \
'Overlap area covers only %s pixels. At least 16*16 pixels are needed.' % px_covered
def equalize_pixGrids(self):
"""
Equalize image grids and projections of reference and target image (align target to reference).
"""
if not (prj_equal(self.ref.prj, self.shift.prj) and self.ref.xygrid_specs==self.shift.xygrid_specs):
if not self.q: print("Equalizing pixel grids and projections of reference and target image...")
if not (prj_equal(self.ref.prj, self.shift.prj) and self.ref.xygrid_specs == self.shift.xygrid_specs):
if not self.q:
print("Equalizing pixel grids and projections of reference and target image...")
if self.grid2use=='ref':
if self.grid2use == 'ref':
# resample target image to refernce image
self.shift.arr = self.shift[:,:,self.shift.band4match] # resample the needed band only
self.shift.arr = self.shift[:, :, self.shift.band4match] # resample the needed band only
self.shift.reproject_to_new_grid(prototype=self.ref, CPUs=self.CPUs)
self.shift.band4match = 0 # after resampling there is only one band in the GeoArray
self.shift.band4match = 0 # after resampling there is only one band in the GeoArray
else:
# resample reference image to target image
# FIXME in case of different projections this will change the projection of the reference image!
self.ref.arr = self.ref[:,:,self.ref.band4match] # resample the needed band only
self.ref.arr = self.ref[:, :, self.ref.band4match] # resample the needed band only
self.ref.reproject_to_new_grid(prototype=self.shift, CPUs=self.CPUs)
self.ref.band4match = 0 # after resampling there is only one band in the GeoArray
def show_image_footprints(self):
"""This method is intended to be called from Jupyter Notebook and shows a web map containing the calculated
footprints of the input images as well as the corresponding overlap area."""
......@@ -448,9 +457,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 .epsg, 4326)
shiftPoly = reproject_shapelyGeometry(self.shift.poly , self.shift.epsg, 4326)
overlapPoly = reproject_shapelyGeometry(self.overlap_poly , self.shift.epsg, 4326)
refPoly = reproject_shapelyGeometry(self.ref.poly, self.ref.epsg, 4326)
shiftPoly = reproject_shapelyGeometry(self.shift.poly, self.shift.epsg, 4326)
overlapPoly = reproject_shapelyGeometry(self.overlap_poly, self.shift.epsg, 4326)
matchBoxPoly = reproject_shapelyGeometry(self.matchBox.mapPoly, self.shift.epsg, 4326)
m = folium.Map(location=tuple(np.array(overlapPoly.centroid.coords.xy).flatten())[::-1])
......@@ -459,8 +468,7 @@ class COREG(object):
folium.GeoJson(gjs).add_to(m)
return m
def show_matchWin(self, figsize=(15,15), interactive=True, after_correction=False):
def show_matchWin(self, figsize=(15, 15), interactive=True, after_correction=False):
"""Show the image content within the matching window.
:param figsize: <tuple> figure size
......@@ -473,54 +481,53 @@ class COREG(object):
try:
import holoviews as hv
except ImportError:
hv =None
hv = None
if not hv:
raise ImportError(
"This method requires the library 'holoviews'. It can be installed for Anaconda with "
"the shell command 'conda install -c ioam holoviews bokeh'.")
warnings.filterwarnings('ignore')
hv.notebook_extension('matplotlib')
hv.Store.add_style_opts(hv.Image, ['vmin','vmax'])
#hv.Store.option_setters.options().Image = hv.Options('style', cmap='gnuplot2')
#hv.Store.add_style_opts(hv.Image, ['cmap'])
#renderer = hv.Store.renderers['matplotlib'].instance(fig='svg', holomap='gif')
#RasterPlot = renderer.plotting_class(hv.Image)
#RasterPlot.cmap = 'gray'
otherWin_corr = self._get_deshifted_otherWin()
xmin,xmax,ymin,ymax = self.matchBox.boundsMap
get_vmin = lambda arr: np.percentile(arr, 2)
get_vmax = lambda arr: np.percentile(arr, 98)
rescale = lambda arr: rescale_intensity(arr, in_range=(get_vmin(arr), get_vmax(arr)))
get_arr = lambda geoArr: rescale(np.ma.masked_equal(geoArr[:], geoArr.nodata))
get_hv_image = lambda geoArr: hv.Image(get_arr(geoArr), bounds=(xmin,ymin,xmax,ymax))(
style={'cmap':'gray',
'vmin':get_vmin(geoArr[:]), 'vmax':get_vmax(geoArr[:]), # does not work
'interpolation':'none'},
plot={'fig_inches':figsize, 'show_grid':True})
#plot={'fig_size':100, 'show_grid':True})
imgs_orig = {1 : get_hv_image(self.matchWin), 2 : get_hv_image(self.otherWin)}
imgs_corr = {1 : get_hv_image(self.matchWin), 2 : get_hv_image(otherWin_corr)}
#layout = get_hv_image(self.matchWin) + get_hv_image(self.otherWin)
imgs = {1 : get_hv_image(self.matchWin) + get_hv_image(self.matchWin),
2 : get_hv_image(self.otherWin) + get_hv_image(otherWin_corr)
}
hv.Store.add_style_opts(hv.Image, ['vmin', 'vmax'])
# hv.Store.option_setters.options().Image = hv.Options('style', cmap='gnuplot2')
# hv.Store.add_style_opts(hv.Image, ['cmap'])
# renderer = hv.Store.renderers['matplotlib'].instance(fig='svg', holomap='gif')
# RasterPlot = renderer.plotting_class(hv.Image)
# RasterPlot.cmap = 'gray'
otherWin_corr = self._get_deshifted_otherWin()
xmin, xmax, ymin, ymax = self.matchBox.boundsMap
get_vmin = lambda arr: np.percentile(arr, 2)
get_vmax = lambda arr: np.percentile(arr, 98)
rescale = lambda arr: rescale_intensity(arr, in_range=(get_vmin(arr), get_vmax(arr)))
get_arr = lambda geoArr: rescale(np.ma.masked_equal(geoArr[:], geoArr.nodata))
get_hv_image = lambda geoArr: hv.Image(get_arr(geoArr), bounds=(xmin, ymin, xmax, ymax))(
style={'cmap': 'gray',
'vmin': get_vmin(geoArr[:]), 'vmax': get_vmax(geoArr[:]), # does not work
'interpolation': 'none'},
plot={'fig_inches': figsize, 'show_grid': True})
# plot={'fig_size':100, 'show_grid':True})
imgs_orig = {1: get_hv_image(self.matchWin), 2: get_hv_image(self.otherWin)}
imgs_corr = {1: get_hv_image(self.matchWin), 2: get_hv_image(otherWin_corr)}
# layout = get_hv_image(self.matchWin) + get_hv_image(self.otherWin)
imgs = {1: get_hv_image(self.matchWin) + get_hv_image(self.matchWin),
2: get_hv_image(self.otherWin) + get_hv_image(otherWin_corr)
}
# Construct a HoloMap by evaluating the function over all the keys
hmap_orig = hv.HoloMap(imgs_orig, kdims=['image'])
hmap_corr = hv.HoloMap(imgs_corr, kdims=['image'])
hmap = hv.HoloMap(imgs, kdims=['image']).collate().cols(1) # displaying this results in a too small figure
#hmap = hv.HoloMap(imgs_corr, kdims=['image']) + hv.HoloMap(imgs_corr, kdims=['image'])
hmap = hv.HoloMap(imgs, kdims=['image']).collate().cols(1) # displaying this results in a too small figure
# hmap = hv.HoloMap(imgs_corr, kdims=['image']) + hv.HoloMap(imgs_corr, kdims=['image'])
## Construct a HoloMap by defining the sampling on the Dimension
#dmap = hv.DynamicMap(image_slice, kdims=[hv.Dimension('z_axis', values=keys)])
# Construct a HoloMap by defining the sampling on the Dimension
# dmap = hv.DynamicMap(image_slice, kdims=[hv.Dimension('z_axis', values=keys)])
warnings.filterwarnings('default')
#return hmap
# return hmap