Commit 216c2693 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

implemented new module 'CoReg_local' for detection and correction of local...

implemented new module 'CoReg_local' for detection and correction of local shifts; revised some functions due to changes in external library 'py_tools_ds'

components.CoReg:
- imParamObj: updated derivation of nodata value and footprint polygo
- COREG: - added attribute 'deshift_results'

components.CoReg_local:
- added as new module for detection and correction of local shifts

components.DeShifter:
- moved DESHIFTER property '_dict_rspAlg_rsp_Int' to module top-level
- _get_out_extent(): updated algorithm

components. Geom_Quality_Grid.Geom_Quality_Grid:
- __init__():
    - changed most of the input args and kwargs
    - revised the whole function -> COREG_obj is not instanced within class COREG_LOCAL
    - added property 'CoRegPoints_table'
    - added property 'GCPList'
    - dump_CoRegPoints_table(): revision
    - quality_grid_to_PointShapefile(): refactored to 'to_PointShapefile'
    - _quality_grid_to_PointShapefile(): refactored to '_to_PointShapefile'
    - quality_grid_to_Raster_using_KrigingOLD(): revised getter for path_out
    - _Kriging_sp(): revised getter for path_out
    - view_results(): moved to class COREG_LOCAL
    - view_results_folium(): moved to class COREG_LOCAL
    - correct_shifts(): moved to class COREG_LOCAL

components.geometry:
- get_true_corner_mapXY(): adjusted some lines due to update of 'py_tools_ds'
- find_noDataVal(): moved to external library 'py_tools_ds'

__init__:
- updated __version__
parent d017de6f
......@@ -8,7 +8,7 @@ from .components import utilities
from .components import geometry
__author__ = 'Daniel Scheffler'
__version__= '2016-10-29_01'
__version__= '2016-11-01_01'
__all__=['COREG',
'DESHIFTER',
......
......@@ -10,7 +10,10 @@ import warnings
from copy import copy
# custom
import gdal
try:
import gdal
except ImportError:
from osgeo import gdal
import numpy as np
try:
import pyfftw
......@@ -19,15 +22,14 @@ except ImportError:
from shapely.geometry import Point
# internal modules
from .DeShifter import DESHIFTER
from .DeShifter import DESHIFTER, _dict_rspAlg_rsp_Int
from . import geometry as GEO
from . import io as IO
from . import plotting as PLT
from py_tools_ds.ptds import GeoArray
from py_tools_ds.ptds.geo.coord_calc import corner_coord_to_minmax, get_corner_coordinates
from py_tools_ds.ptds.geo.vector.topology import get_footprint_polygon, get_overlap_polygon, \
get_smallest_boxImYX_that_contains_boxMapYX
from py_tools_ds.ptds.geo.vector.topology import get_overlap_polygon, get_smallest_boxImYX_that_contains_boxMapYX
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
......@@ -79,7 +81,7 @@ class imParamObj(object):
if CoReg_params['nodata'][0 if imID == 'ref' else 1] is not None:
self.nodata = CoReg_params['nodata'][0 if imID == 'ref' else 1]
else:
self.nodata = GEO.find_noDataVal(self.GeoArray)
self.nodata = self.GeoArray.nodata
# set corner coords
given_corner_coord = CoReg_params['data_corners_%s' % ('im0' if imID == 'ref' else 'im1')]
......@@ -96,7 +98,8 @@ class imParamObj(object):
self.corner_coord = given_corner_coord
# set footprint polygon
self.poly = get_footprint_polygon(self.corner_coord, fix_invalid=True)
#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)), \
"The corner position '%s' is outside of the %s." % (XY, self.imName)
......@@ -105,8 +108,6 @@ class imParamObj(object):
class COREG(object):
_dict_rspAlg_rsp_Int = DESHIFTER._dict_rspAlg_rsp_Int
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,
......@@ -171,7 +172,7 @@ class COREG(object):
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 self._dict_rspAlg_rsp_Int.keys(), "'%s' is not a supported resampling algorithm." % rspAlg
assert rspAlg in _dict_rspAlg_rsp_Int.keys(), "'%s' is not a supported resampling algorithm." % rspAlg
if resamp_alg_calc=='average':
warnings.warn("The resampling algorithm 'average' causes sinus-shaped patterns in fft images that will "
"affect the precision of the calculated spatial shifts! It is highly recommended to"
......@@ -214,6 +215,7 @@ class COREG(object):
self.tracked_errors = [] # expanded each time an error occurs
self.success = False # default
self.deshift_results = None # set by self.correct_shifts()
gdal.AllRegister()
self._get_image_params()
......@@ -360,8 +362,18 @@ class COREG(object):
overlap_center_pos_x, overlap_center_pos_y = self.overlap_poly.centroid.coords.xy
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.'
#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
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)
......@@ -848,8 +860,8 @@ class COREG(object):
CPUs = None if self.mp else 1,
v = self.v,
q = self.q)
deshift_results = DS.correct_shifts()
return deshift_results
self.deshift_results = DS.correct_shifts()
return self.deshift_results
def correct_shifts_OLD(self):
......
# -*- coding: utf-8 -*-
__author__='Daniel Scheffler'
import warnings
import os
# custom
try:
import gdal
except ImportError:
from osgeo import gdal
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from .Geom_Quality_Grid import Geom_Quality_Grid
from .CoReg import COREG
from .DeShifter import DESHIFTER
from py_tools_ds.ptds.geo.coord_trafo import transform_any_prj, reproject_shapelyGeometry
from py_tools_ds.ptds import GeoArray
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):
"""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
correction performs a polynomial transformation using te calculated shifts of each point in the grid as GCPs.
Thus 'Geom_Quality_Grid' can be used to correct for locally varying geometric distortions of the target image.
:param im_ref(str, GeoArray): source path of reference image (any GDAL compatible image format is supported)
:param im_tgt(str, GeoArray): source path of image to be shifted (any GDAL compatible image format is supported)
:param grid_res: grid resolution in pixels of the target image
:param window_size(tuple): custom matching window size [pixels] (default: (512,512))
:param path_out(str): target path of the coregistered image
- if None (default), the method correct_shifts() does not write to disk
- if 'auto': /dir/of/im1/<im1>__shifted_to__<im0>.bsq
:param fmt_out(str): raster file format for output file. ignored if path_out is None. can be any GDAL
compatible raster file format (e.g. 'ENVI', 'GeoTIFF'; default: ENVI)
:param projectDir:
:param r_b4match(int): band of reference image to be used for matching (starts with 1; default: 1)
: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 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
: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 binary_ws(bool): use binary X/Y dimensions for the matching window (default: 1)
: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)
"""
self.params = dict([x for x in locals().items() if x[0] != "self" and not x[0].startswith('__')])
self.imref = im_ref if isinstance(im_ref, GeoArray) else GeoArray(im_ref)
self.im2shift = im_tgt if isinstance(im_tgt, GeoArray) else GeoArray(im_tgt)
self.path_out = path_out # updated by self.set_outpathes
self.fmt_out = fmt_out
self._projectDir = projectDir
self.grid_res = grid_res
self.window_size = window_size
self.max_shift = max_shift
self.max_iter = max_iter
self.r_b4match = r_b4match
self.s_b4match = s_b4match
self.calc_corners = calc_corners
self.nodata = nodata
self.outFillVal = outFillVal
self.bin_ws = binary_ws
self.CPUs = CPUs
self.v = v
self.q = q
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
if path_out and projectDir and os.path.basename(self.path_out):
self.path_out = os.path.join(self.projectDir, os.path.basename(self.path_out))
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,
v = v,
q = q,
ignore_errors = True)
self._quality_grid = None # set by self.quality_grid
self._CoRegPoints_table = None # set by self.CoRegPoints_table
self.deshift_results = None # set by self.correct_shifts()
@property
def projectDir(self):
if self._projectDir:
if len(os.path.split(self._projectDir))==1:
return os.path.abspath(os.path.join(os.path.curdir, self._projectDir))
else:
return os.path.abspath(self._projectDir)
else:
# return a project name that not already has a corresponding folder on disk
projectDir = os.path.join(os.path.dirname(self.im2shift.filePath), 'UntitledProject_1')
while os.path.isdir(projectDir):
projectDir = '%s_%s' % (projectDir.split('_')[0], int(projectDir.split('_')[1]) + 1)
self._projectDir = projectDir
return self._projectDir
@property
def quality_grid(self):
if self._quality_grid:
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)
return self._quality_grid
@property
def CoRegPoints_table(self):
return self.quality_grid.CoRegPoints_table
def view_CoRegPoints(self, attribute2plot='ABS_SHIFT', cmap=None, exclude_fillVals=True, backgroundIm='tgt',
figsize=None, savefigPath='', savefigDPI=96):
"""Shows a map of the calculated quality grid with the target image as background.
:param attribute2plot: <str> the attribute of the quality grid to be shown (default: 'ABS_SHIFT')
:param cmap: <plt.cm.<colormap>> a custom color map to be applied to the plotted grid points
(default: 'RdYlGn_r')
:param exclude_fillVals: <bool> whether to exclude those points of the grid where spatial shift detection failed
:param backgroundIm: <str> whether to use the target or the reference image as map background. Possible
options are 'ref' and 'tgt' (default: 'tgt')
:param figsize: <tuple> size of the figure to be viewed, e.g. (10,10)
:param savefigPath:
:param savefigDPI:
:return:
"""
# get a map showing target image
if backgroundIm not in ['tgt','ref']: raise ValueError('backgroundIm')
backgroundIm = self.im2shift if backgroundIm=='tgt' else self.imref
fig, ax, map2show = backgroundIm.show_map(figsize=figsize, nodataVal=self.nodata[1], return_map=True)
# fig, ax, map2show = backgroundIm.show_map_utm(figsize=(20,20), nodataVal=self.nodata[1], return_map=True)
plt.title(attribute2plot)
# transform all points of quality grid to LonLat
GDF = self.CoRegPoints_table.loc[self.CoRegPoints_table.X_SHIFT_M != self.outFillVal, ['geometry', attribute2plot]].copy() \
if exclude_fillVals else self.CoRegPoints_table.loc[:, ['geometry', attribute2plot]]
# get LonLat coordinates for all points
get_LonLat = lambda X, Y: transform_any_prj(self.im2shift.projection, 4326, X, Y)
GDF['LonLat'] = list(GDF['geometry'].map(lambda geom: get_LonLat(*tuple(np.array(geom.coords.xy)[:,0]))))
# get colors for all points
#vmin = min(GDF[GDF[attribute2plot] != self.outFillVal][attribute2plot])
#vmax = max(GDF[GDF[attribute2plot] != self.outFillVal][attribute2plot])
#norm = mpl_normalize(vmin=vmin, vmax=vmax)
palette = cmap if cmap else plt.cm.RdYlGn_r
#GDF['color'] = [*GDF[attribute2plot].map(lambda val: palette(norm(val)))]
# add quality grid to map
#plot_point = lambda row: ax.plot(*map2show(*row['LonLat']), marker='o', markersize=7.0, alpha=1.0, color=row['color'])
#GDF.apply(plot_point, axis=1)
GDF['plt_XY'] = list(GDF['LonLat'].map(lambda ll: map2show(*ll)))
GDF['plt_X'] = list(GDF['plt_XY'].map(lambda XY: XY[0]))
GDF['plt_Y'] = list(GDF['plt_XY'].map(lambda XY: XY[1]))
points = plt.scatter(GDF['plt_X'],GDF['plt_Y'], c=GDF[attribute2plot],
cmap=palette, marker='o' if len(GDF)<10000 else '.', s=50, alpha=1.0)
# add colorbar
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", size="2%", pad=0.1) # create axis on the right; size =2% of ax; padding = 0.1 inch
plt.colorbar(points, cax=cax)
plt.show()
if savefigPath:
fig.savefig(savefigPath, dpi=savefigDPI)
def view_CoRegPoints_folium(self, attribute2plot='ABS_SHIFT', cmap=None, exclude_fillVals=True):
warnings.warn(UserWarning('This function is still under construction and may not work as expected!'))
assert self.CoRegPoints_table is not None, 'Calculate quality grid first!'
try:
import folium, geojson
from folium import plugins
except ImportError:
folium, geojson, plugins = [None]*3
if not folium or not geojson:
raise ImportError("This method requires the libraries 'folium' and 'geojson'. They can be installed with "
"the shell command 'pip install folium geojson'.")
lon_min, lat_min, lon_max, lat_max = \
reproject_shapelyGeometry(self.im2shift.box.mapPoly, self.im2shift.projection, 4326).bounds
center_lon, center_lat = (lon_min+lon_max)/2, (lat_min+lat_max)/2
# get image to plot
image2plot = self.im2shift[:]
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],
out_XYdims=(1000, 1000), q=True, out_prj='epsg:3857') # image must be transformed into web mercator projection
# create map
map_osm = folium.Map(location=[center_lat, center_lon])#,zoom_start=3)
import matplotlib
plugins.ImageOverlay(
colormap=lambda x: (1, 0, 0, x), # TODO a colormap must be given
# colormap=matplotlib.cm.gray, # does not work
image=image2plot, bounds=[[lat_min, lon_min], [lat_max, lon_max]],
).add_to(map_osm)
folium.GeoJson(self.CoRegPoints_table.loc[:, ['geometry', attribute2plot]]).add_to(map_osm)
# add overlap polygon
overlapPoly = reproject_shapelyGeometry(self.COREG_obj.overlap_poly, self.im2shift.epsg, 4326)
gjs = geojson.Feature(geometry=overlapPoly, properties={})
folium.GeoJson(gjs).add_to(map_osm)
return map_osm
def correct_shifts(self, max_GCP_count=None):
"""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
:return:
"""
coreg_info = self.COREG_obj.coreg_info
coreg_info['GCPList'] = self.quality_grid.GCPList
if max_GCP_count:
coreg_info['GCPList'] = coreg_info['GCPList'][:max_GCP_count] # TODO should be a random sample
DS = DESHIFTER(self.im2shift, coreg_info,
path_out = self.path_out,
fmt_out = self.fmt_out,
out_gsd = (self.im2shift.xgsd,self.im2shift.ygsd),
align_grids = True,
#cliptoextent = True, # why?
#cliptoextent = False,
#clipextent = self.im2shift.box.boxMapYX,
#options = '-wm 10000 -order 1',
v = self.v,
q = self.q)
self.deshift_results = DS.correct_shifts()
return self.deshift_results
......@@ -8,24 +8,24 @@ import time
import warnings
# custom
import gdal
try:
import gdal
except ImportError:
from osgeo import gdal
import numpy as np
import rasterio
from shapely.geometry import box
# internal modules
from . import geometry as GEO
from py_tools_ds.ptds import GeoArray
from py_tools_ds.ptds.numeric.array import get_outFillZeroSaturated
from py_tools_ds.ptds.geo.map_info import mapinfo2geotransform, geotransform2mapinfo
from py_tools_ds.ptds.geo.coord_calc import corner_coord_to_minmax, get_corner_coordinates
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}
class DESHIFTER(object):
def __init__(self, im2shift, coreg_results, **kwargs):
......@@ -82,7 +82,7 @@ class DESHIFTER(object):
self.path_out = kwargs.get('path_out' , None)
self.fmt_out = kwargs.get('fmt_out' , 'ENVI')
self.band2process = kwargs.get('band2process', None) # starts with 1 # FIXME warum?
self.nodata = kwargs.get('nodata' , get_outFillZeroSaturated(self.im2shift.dtype)[0]) # TODO auto-detection here?
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
......@@ -98,7 +98,7 @@ class DESHIFTER(object):
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 self._dict_rspAlg_rsp_Int.keys(), \
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']
......@@ -110,12 +110,6 @@ class DESHIFTER(object):
self.GeoArray_shifted = None # set by self.correct_shifts
@property
def _dict_rspAlg_rsp_Int(self):
return {'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}
def _get_out_grid(self, init_kwargs):
# parse given params
out_gsd = init_kwargs.get('out_gsd' , None)
......@@ -181,13 +175,9 @@ class DESHIFTER(object):
def _get_out_extent(self):
if self.cliptoextent and self.clipextent is None:
# calculate actual corner coords (input image projection)
trueCorner_mapXY = GEO.get_true_corner_mapXY(self.im2shift, bandNr=1, noDataVal=self.nodata, mp=1,v=0)
xmin, xmax, ymin, ymax = corner_coord_to_minmax(trueCorner_mapXY)
self.clipextent = box(xmin, ymin, xmax, ymax).bounds
self.clipextent = self.im2shift.footprint_poly.bounds
else:
xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(
gt=self.shift_gt, cols=self.im2shift.cols, rows=self.im2shift.rows))
xmin, xmax, ymin, ymax = self.im2shift.box.boundsMap
self.clipextent = xmin, ymin, xmax, ymax
......@@ -285,13 +275,13 @@ class DESHIFTER(object):
# get resampled array
out_arr, out_gt, out_prj = \
warp_ndarray(in_arr, self.shift_gt, self.shift_prj, self.ref_prj,
rspAlg = self._dict_rspAlg_rsp_Int[self.rspAlg],
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,
polynomialOrder = None,
options = None, #'-refine_gcps 500',
CPUs = self.CPUs,
q = self.q)
......
This diff is collapsed.
......@@ -18,7 +18,6 @@ except ImportError:
from osgeo import ogr
# internal modules
from py_tools_ds.ptds import GeoArray
from py_tools_ds.ptds.geo.coord_calc import calc_FullDataset_corner_positions
from py_tools_ds.ptds.geo.coord_trafo import pixelToMapYX, imYX2mapYX
from py_tools_ds.ptds import GeoArray
......@@ -41,15 +40,16 @@ def get_true_corner_mapXY(fPath_or_geoarray, bandNr=1, noDataVal=None, mp=1, v=0
rows,cols = geoArr.shape[:2]
gt, prj = geoArr.geotransform, geoArr.projection
assert gt and prj, 'GeoTransform an projection must be given for calculation of LonLat corner coordinates.'
mask_1bit = np.zeros((rows,cols),dtype='uint8') # zeros -> image area later overwritten by ones
if noDataVal is None:
mask_1bit[:,:] = 1
elif noDataVal=='unclear':
elif noDataVal=='ambiguous':
warnings.warn("No data value could not be automatically detected. Thus the matching window used for shift "
"calulation had to be centered in the middle of the overlap center without respecting no data values. "
"calculation had to be centered in the middle of the overlap area without respecting no data values. "
"To avoid this provide the correct no data values for reference and shift image via '-nodata'")
mask_1bit[:,:] = 1
else:
......@@ -64,7 +64,7 @@ def get_true_corner_mapXY(fPath_or_geoarray, bandNr=1, noDataVal=None, mp=1, v=0
if v:
warnings.warn("\nCalculation of corner coordinates failed within algorithm 'shapely' (Exception: %s)."
" Using algorithm 'numpy' instead." %sys.exc_info()[1])
corner_coords_YX = calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=False, algorithm='numpy')
corner_coords_YX = calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=False, algorithm='numpy') # FIXME numpy algorithm returns wrong values for S2A_OPER_MSI_L1C_TL_SGS__20160608T153121_A005024_T33UUU_B03.jp2 (Hannes)
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
......@@ -80,7 +80,7 @@ def get_true_corner_mapXY(fPath_or_geoarray, bandNr=1, noDataVal=None, mp=1, v=0
#all_coords_are_unique = len([UL, UR, LL, LR]) == len(GeoDataFrame([UL, UR, LL, LR]).drop_duplicates().values)
#UL, UR, LL, LR = (UL, UR, LL, LR) if all_coords_are_unique else ((0, 0), (0, cols-1), (rows-1, 0), (rows-1, cols-1))
get_mapYX = lambda YX: pixelToMapYX(list(reversed(YX)),geotransform=geoArr.geotransform,projection=prj)[0]
get_mapYX = lambda YX: pixelToMapYX(list(reversed(YX)),geotransform=gt,projection=prj)[0]
corner_pos_XY = [list(reversed(i)) for i in [get_mapYX(YX) for YX in corner_coords_YX]]
return corner_pos_XY
......@@ -103,21 +103,4 @@ def get_GeoArrayPosition_from_boxImYX(boxImYX):
"""Returns row_start,row_end,col_start,col_end and assumes boxImYX as [UL_YX,UR_YX,LR_YX,LL_YX)"""
rS, cS = boxImYX[0] # UL
rE, cE = boxImYX[2] # LR
return rS, rE, cS, cE
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 pathIm_or_GeoArray:
:param bandIdx:
:param sz: window size in which corner pixels are analysed
"""
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)}
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
return rS, rE, cS, cE
\ 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