Commit 82509a5a authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

some further developments and bug fixes

io.raster.gdal:
- catched empty dataset exception

io.raster.GeoArray.GeoArray:
- added property 'xygrid_specs'
- save(): edited docstring; catched missing driver exception
- _get_plottable_image(): now supports requesting a geographic area (added keywords 'boundsMap', 'boundsMapPrj')
- show(): implemented keywords 'boundsMap', 'boundsMapPrj'; edited docstring
- show_map(): implemented keywords 'boundsMap', 'boundsMapPrj'; edited docstring
- get_mapPos(): refactored keyword bandslist to band2get; edited docstring

numeric.vector:
- find_nearest(): bugfix for wrong return value in case of exact hit if roundAlg=='off'

similarity.raster:
- calc_ssim(): added data type adjustment; refactored match to image0 and other to image1; added keyword 'gaussian_weights'

- updated __version__
parent 44ce3f7a
...@@ -15,7 +15,7 @@ __all__=[#'compatibility', ...@@ -15,7 +15,7 @@ __all__=[#'compatibility',
'similarity', 'similarity',
'GeoArray'] 'GeoArray']
__version__ = '20161182_02' __version__ = '20161112_01'
__author__='Daniel Scheffler' __author__='Daniel Scheffler'
# Validate GDAL version # Validate GDAL version
......
...@@ -211,6 +211,14 @@ class GeoArray(object): ...@@ -211,6 +211,14 @@ class GeoArray(object):
return abs(self.geotransform[5]) return abs(self.geotransform[5])
@property
def xygrid_specs(self):
"""Get the specifications for the X/Y coordinate grid, e.g. [[15,30], [0,30]] for a coordinate with its origin
at X/Y[15,0] and a GSD of X/Y[15,30]."""
get_grid = lambda gt, xgsd, ygsd: [[gt[0], gt[0] + xgsd], [gt[3], gt[3] - ygsd]]
return get_grid(self.geotransform, self.xgsd, self.ygsd)
@property @property
def projection(self): def projection(self):
"""Get the projection of the associated image. Setting the projection is only allowed if GeoArray has been """Get the projection of the associated image. Setting the projection is only allowed if GeoArray has been
...@@ -565,11 +573,12 @@ class GeoArray(object): ...@@ -565,11 +573,12 @@ class GeoArray(object):
def save(self, out_path, fmt='ENVI', creationOptions=None): def save(self, out_path, fmt='ENVI', creationOptions=None):
# type: (str, str, list)
"""Write the raster data to disk. """Write the raster data to disk.
:param out_path: <str> output path :param out_path: <str> output path
:param fmt: <str> the output format / GDAL driver code to be used for output creation, e.g. 'ENVI' :param fmt: <str> the output format / GDAL driver code to be used for output creation, e.g. 'ENVI'
:param creationOptions: <list> GDAL creation options, e.g. ["QUALITY=20", "REVERSIBLE=YES", "WRITE_METADATA=YES"] :param creationOptions: <list> GDAL creation options, e.g. ["QUALITY=80", "REVERSIBLE=YES", "WRITE_METADATA=YES"]
:return: :return:
""" """
...@@ -579,8 +588,11 @@ class GeoArray(object): ...@@ -579,8 +588,11 @@ class GeoArray(object):
if not os.path.isdir(os.path.dirname(out_path)): os.makedirs(os.path.dirname(out_path)) if not os.path.isdir(os.path.dirname(out_path)): os.makedirs(os.path.dirname(out_path))
if self.is_inmem: if self.is_inmem:
ds = get_GDAL_ds_inmem(self.arr,self.geotransform, self.projection, self.nodata) # expects rows,columns,bands ds = get_GDAL_ds_inmem(self.arr,self.geotransform, self.projection, self.nodata) # expects rows,columns,bands
gdal.GetDriverByName(fmt).CreateCopy(out_path, ds, options=creationOptions if creationOptions else []) driver = gdal.GetDriverByName(fmt)
if driver is None:
raise Exception("'%s' is not a supported GDAL driver." %fmt)
driver.CreateCopy(out_path, ds, options=creationOptions if creationOptions else [])
#out_arr = self.arr if self.ndim == 2 else np.swapaxes(np.swapaxes(self.arr, 0, 2), 1, 2) # rows, columns, bands => bands, rows, columns #out_arr = self.arr if self.ndim == 2 else np.swapaxes(np.swapaxes(self.arr, 0, 2), 1, 2) # rows, columns, bands => bands, rows, columns
#gdalnumeric.SaveArray(out_arr, out_path, format=fmt, prototype=ds) # expects bands,rows,columns #gdalnumeric.SaveArray(out_arr, out_path, format=fmt, prototype=ds) # expects bands,rows,columns
...@@ -592,21 +604,33 @@ class GeoArray(object): ...@@ -592,21 +604,33 @@ class GeoArray(object):
gdal_Translate(out_path, src_ds, format=fmt, creationOptions=creationOptions) gdal_Translate(out_path, src_ds, format=fmt, creationOptions=creationOptions)
src_ds = None src_ds = None
if not os.path.exists(out_path):
raise Exception(gdal.GetLastErrorMsg())
def dump(self, out_path): def dump(self, out_path):
# type: (str)
"""Sertialize the whole object instance to disk using dill.""" """Sertialize the whole object instance to disk using dill."""
import dill import dill
with open(out_path,'w') as outF: with open(out_path,'w') as outF:
dill.dump(self,outF) dill.dump(self,outF)
def _get_plottable_image(self, xlim=None, ylim=None, band=0, res_factor=None, nodataVal=None, out_prj=None): def _get_plottable_image(self, xlim=None, ylim=None, band=0, boundsMap=None, boundsMapPrj=None, res_factor=None,
nodataVal=None, out_prj=None):
# handle limits # handle limits
cS, cE = xlim if isinstance(xlim, (tuple, list)) else (0, self.columns - 1) if boundsMap:
rS, rE = ylim if isinstance(ylim, (tuple, list)) else (0, self.rows - 1) boundsMapPrj = boundsMapPrj if boundsMapPrj else self.prj
image2plot, gt, prj = self.get_mapPos(boundsMap, boundsMapPrj, band2get=band,
fillVal=nodataVal if nodataVal is not None else self.nodata)
else:
cS, cE = xlim if isinstance(xlim, (tuple, list)) else (0, self.columns - 1)
rS, rE = ylim if isinstance(ylim, (tuple, list)) else (0, self.rows - 1)
image2plot = self[rS:rE, cS:cE, band]
gt, prj = self.geotransform, self.projection
image2plot = self[rS:rE, cS:cE, band]
gt, prj = self.geotransform, self.projection
transOpt = ['SRC_METHOD=NO_GEOTRANSFORM'] if tuple(gt) == (0, 1, 0, 0, 0, -1) else None transOpt = ['SRC_METHOD=NO_GEOTRANSFORM'] if tuple(gt) == (0, 1, 0, 0, 0, -1) else None
xdim, ydim = None, None xdim, ydim = None, None
nodataVal = nodataVal if nodataVal is not None else self.nodata nodataVal = nodataVal if nodataVal is not None else self.nodata
...@@ -633,13 +657,15 @@ class GeoArray(object): ...@@ -633,13 +657,15 @@ class GeoArray(object):
return image2plot, gt, prj return image2plot, gt, prj
def show(self, xlim=None, ylim=None, band=0, figsize=None, interpolation='none', cmap=None, nodataVal=None, def show(self, xlim=None, ylim=None, band=0, boundsMap=None, boundsMapPrj=None, figsize=None, interpolation='none',
res_factor=None): cmap=None, nodataVal=None, res_factor=None):
"""Plots the desired array position into a figure. """Plots the desired array position into a figure.
:param xlim: :param xlim: [start_column, end_column]
:param ylim: :param ylim: [start_row, end_row]
:param band: :param band:
:param boundsMap: xmin, ymin, xmax, ymax
:param boundsMapPrj:
:param figsize: :param figsize:
:param interpolation: :param interpolation:
:param cmap: :param cmap:
...@@ -650,7 +676,9 @@ class GeoArray(object): ...@@ -650,7 +676,9 @@ class GeoArray(object):
# get image to plot # get image to plot
nodataVal = nodataVal if nodataVal is not None else self.nodata nodataVal = nodataVal if nodataVal is not None else self.nodata
image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal) image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, boundsMap=boundsMap,
boundsMapPrj=boundsMapPrj, res_factor=res_factor,
nodataVal=nodataVal)
# set color palette # set color palette
palette = cmap if cmap else plt.cm.gray palette = cmap if cmap else plt.cm.gray
...@@ -665,18 +693,21 @@ class GeoArray(object): ...@@ -665,18 +693,21 @@ class GeoArray(object):
# show image # show image
plt.figure(figsize=figsize) plt.figure(figsize=figsize)
plt.imshow(image2plot, palette, interpolation=interpolation, extent=(0, self.columns, self.rows, 0), rows, cols = image2plot.shape[:2]
plt.imshow(image2plot, palette, interpolation=interpolation, extent=(0, cols, rows, 0),
vmin=vmin, vmax=vmax, ) # compressed excludes nodata values vmin=vmin, vmax=vmax, ) # compressed excludes nodata values
plt.show() plt.show()
def show_map(self, xlim=None, ylim=None, band=0, figsize=None, interpolation='none', cmap=None, nodataVal=None, def show_map(self, xlim=None, ylim=None, band=0, boundsMap=None, boundsMapPrj=None, figsize=None,
res_factor=None, return_map=False): interpolation='none', cmap=None, nodataVal=None, res_factor=None, return_map=False):
""" """
:param xlim: :param xlim:
:param ylim: :param ylim:
:param band: :param band:
:param boundsMap: xmin, ymin, xmax, ymax
:param boundsMapPrj:
:param figsize: :param figsize:
:param interpolation: :param interpolation:
:param cmap: :param cmap:
...@@ -693,9 +724,18 @@ class GeoArray(object): ...@@ -693,9 +724,18 @@ class GeoArray(object):
# get image to plot # get image to plot
nodataVal = nodataVal if nodataVal is not None else self.nodata nodataVal = nodataVal if nodataVal is not None else self.nodata
image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal, 'epsg:4326') image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, boundsMap=boundsMap,
boundsMapPrj=boundsMapPrj, res_factor=res_factor,
nodataVal=nodataVal, out_prj='epsg:4326')
# calculate corner coordinates of plot # calculate corner coordinates of plot
#if boundsMap:
# boundsMapPrj = boundsMapPrj if boundsMapPrj else self.prj
# if not prj_equal(boundsMapPrj, 4326):
# boundsMap = reproject_shapelyGeometry(box(*boundsMap), boundsMapPrj, 4626).bounds
# xmin, ymin, xmax, ymax = boundsMap
# UL_XY, UR_XY, LR_XY, LL_XY = (xmin,ymax), (xmax, ymax), (xmax,ymin), (xmin, ymin)
#else:
UL_XY, UR_XY, LR_XY, LL_XY = [(YX[1],YX[0]) for YX in GeoArray(image2plot, gt, prj).box.boxMapYX] UL_XY, UR_XY, LR_XY, LL_XY = [(YX[1],YX[0]) for YX in GeoArray(image2plot, gt, prj).box.boxMapYX]
center_lon, center_lat = (UL_XY[0]+UR_XY[0])/2., (UL_XY[1]+LL_XY[1])/2. center_lon, center_lat = (UL_XY[0]+UR_XY[0])/2., (UL_XY[1]+LL_XY[1])/2.
...@@ -814,15 +854,16 @@ class GeoArray(object): ...@@ -814,15 +854,16 @@ class GeoArray(object):
return m return m
def get_mapPos(self, mapBounds, mapBounds_prj, bandslist=None, arr_gt=None, arr_prj=None, fillVal=0, rspAlg='near'): # TODO implement slice for indexing bands def get_mapPos(self, mapBounds, mapBounds_prj, band2get=None, arr_gt=None, arr_prj=None, fillVal=None,
rspAlg='near'): # TODO implement slice for indexing bands
""" """
:param mapBounds: xmin, ymin, xmax, ymax :param mapBounds: xmin, ymin, xmax, ymax
:param mapBounds_prj: :param mapBounds_prj:
:param bandslist: :param band2get: band index of the band to be returned (full array if not given)
:param arr_gt: :param arr_gt:
:param arr_prj: :param arr_prj:
:param fillVal: :param fillVal: nodata value
:param rspAlg: <str> Resampling method to use. Available methods are: :param rspAlg: <str> Resampling method to use. Available methods are:
near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, q1, q2 near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, q1, q2
:return: :return:
...@@ -830,12 +871,13 @@ class GeoArray(object): ...@@ -830,12 +871,13 @@ class GeoArray(object):
arr_gt = arr_gt if arr_gt else self.geotransform arr_gt = arr_gt if arr_gt else self.geotransform
arr_prj = arr_prj if arr_prj else self.projection arr_prj = arr_prj if arr_prj else self.projection
fillVal = fillVal if fillVal is not None else self.nodata
if self.is_inmem and not arr_gt or not arr_prj: if self.is_inmem and not arr_gt or not arr_prj:
raise ValueError('In case of in-mem arrays the respective geotransform and projection of the array ' raise ValueError('In case of in-mem arrays the respective geotransform and projection of the array '
'has to be passed.') 'has to be passed.')
sub_arr, sub_gt, sub_prj = get_array_at_mapPos(self, arr_gt, arr_prj, mapBounds_prj, mapBounds, fillVal=fillVal, sub_arr, sub_gt, sub_prj = get_array_at_mapPos(self, arr_gt, arr_prj, mapBounds_prj, mapBounds, fillVal=fillVal,
rspAlg=rspAlg, out_gsd=(self.xgsd,self.ygsd)) rspAlg=rspAlg, out_gsd=(self.xgsd,self.ygsd), band2get=band2get)
return sub_arr, sub_gt, sub_prj return sub_arr, sub_gt, sub_prj
...@@ -1046,7 +1088,7 @@ def get_array_at_mapPos(arr, arr_gt, arr_prj, out_prj, mapBounds, mapBounds_prj= ...@@ -1046,7 +1088,7 @@ def get_array_at_mapPos(arr, arr_gt, arr_prj, out_prj, mapBounds, mapBounds_prj=
# check if reprojection is needed # check if reprojection is needed
mapBounds_prj = mapBounds_prj if mapBounds_prj else out_prj mapBounds_prj = mapBounds_prj if mapBounds_prj else out_prj
samePrj = prj_equal(arr_prj, out_prj) samePrj = prj_equal(arr_prj, out_prj)
if samePrj: if samePrj:
out_prj = arr_prj out_prj = arr_prj
......
...@@ -30,6 +30,8 @@ def get_GDAL_ds_inmem(array, gt=None, prj=None, nodata=None): ...@@ -30,6 +30,8 @@ def get_GDAL_ds_inmem(array, gt=None, prj=None, nodata=None):
OpenNumPyArray = get_gdalnumeric_func('OpenNumPyArray') OpenNumPyArray = get_gdalnumeric_func('OpenNumPyArray')
ds = OpenNumPyArray(array) ds = OpenNumPyArray(array)
if ds is None:
raise Exception(gdal.GetLastErrorMsg())
if gt: ds.SetGeoTransform(gt) if gt: ds.SetGeoTransform(gt)
if prj: ds.SetProjection(prj) if prj: ds.SetProjection(prj)
if nodata: if nodata:
......
...@@ -41,7 +41,11 @@ def find_nearest(array, value, roundAlg='auto', extrapolate=False, exclude_val=F ...@@ -41,7 +41,11 @@ def find_nearest(array, value, roundAlg='auto', extrapolate=False, exclude_val=F
isMiddleVal = collections.Counter(diffs)[minDiff] > 1 isMiddleVal = collections.Counter(diffs)[minDiff] > 1
out = array[minIdx] if not isMiddleVal else array[bisect.bisect_left(array, value)] out = array[minIdx] if not isMiddleVal else array[bisect.bisect_left(array, value)]
elif roundAlg == 'off': elif roundAlg == 'off':
out = array[bisect.bisect_left(array, value) - 1] idx = bisect.bisect_left(array, value)
if array[idx]==value:
out = value # exact hit
else:
out = array[idx-1] # round off
else: # roundAlg == 'on' else: # roundAlg == 'on'
out = array[bisect.bisect_left(array, value)] out = array[bisect.bisect_left(array, value)]
......
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler" __author__ = "Daniel Scheffler"
import numpy as np
from skimage.measure import compare_ssim as ssim from skimage.measure import compare_ssim as ssim
def calc_ssim(match, other, dynamic_range=None, win_size=None): def calc_ssim(image0, image1, dynamic_range=None, win_size=None, gaussian_weights=False):
"""Calculates Mean Structural Similarity Index between two images. """Calculates Mean Structural Similarity Index between two images.
:param match: :param image0:
:param other: :param image1:
:param dynamic_range: :param dynamic_range:
:param win_size: :param win_size:
:param gaussian_weights:
:return: :return:
""" """
return ssim(match, other, dynamic_range=dynamic_range, win_size=win_size) if image0.dtype != image1.dtype:
\ No newline at end of file image0 = image0.astype(np.int16)
image1 = image1.astype(np.int16)
return ssim(image0, image1, dynamic_range=dynamic_range, win_size=win_size, gaussian_weights=gaussian_weights)
\ No newline at end of file
Supports Markdown
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