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

added function for quickly reprojecting GeoArray to a given pixel grid;...

added function for quickly reprojecting GeoArray to a given pixel grid; mask_nodata is now a subclass of GeoArray

geo.coord_grid:
- snap_bounds_to_pixGrid(): added docstring

io.raster.GeoArray:
-GeoArray:
    - bandnames getter: bugfix for not checking if shape of GeoArray.arr changed
    - mask_nodata setter: now sets mask_nodata to a subclass of GeoArray 'NoDataMask'
    - __getattr__(): now also returns results of functions belonging to np.array instances
    - calc_mask_nodata(): changed data type of nodata mask to np.bool
    - save(): bugfix for wrong metadata keys
    - added reproject_to_new_grid(): function for quickly reprojection all array-like attributes to a given target grid
- BadDataMask:
    - bugfix for not allowing certain pixel value combinations
added class NoDataMask

updated __version__
parent 7896b92e
......@@ -15,7 +15,7 @@ __all__=[#'compatibility',
'similarity',
'GeoArray']
__version__ = '20161122_01'
__version__ = '20161125_01'
__author__='Daniel Scheffler'
# Validate GDAL version
......
......@@ -10,6 +10,14 @@ from .coord_calc import get_corner_coordinates
def snap_bounds_to_pixGrid(bounds, gt, roundAlg='auto'):
"""Snaps the given bounds to the given grid (defined by gt) under the use of the given round algorithm.
NOTE: asserts equal projections of source and target grid
:param bounds: <tuple, list> (xmin, ymin, xmax, ymax)
:param gt: <tuple, list> GDAL geotransform
:param roundAlg: <str> 'auto', 'on', 'off'
:return:
"""
in_xmin, in_ymin, in_xmax, in_ymax = bounds
xgrid = np.arange(gt[0], gt[0] + 2 * gt[1], gt[1])
ygrid = np.arange(gt[3], gt[3] + 2 * abs(gt[5]), abs(gt[5]))
......
......@@ -138,7 +138,7 @@ class GeoArray(object):
@property
def bandnames(self):
if self._bandnames:
if self._bandnames and len(self._bandnames)==self.bands:
return self._bandnames
else:
self._bandnames = {'B%s' % band: i for i, band in enumerate(range(1, self.bands + 1))}
......@@ -340,8 +340,30 @@ class GeoArray(object):
@mask_nodata.setter
def mask_nodata(self, maskarray):
self._mask_nodata = maskarray.astype(np.uint8)
def mask_nodata(self, mask):
"""Set bad data mask.
:param mask: Can be a file path, a numpy array or an instance o GeoArray.
"""
if mask is not None:
geoArr_mask = NoDataMask(mask)
geoArr_mask.gt = geoArr_mask.gt if geoArr_mask.gt not in [None, [0, 1, 0, 0, 0, -1]] else self.gt
geoArr_mask.prj = geoArr_mask.prj if geoArr_mask.prj is None else self.prj
imName = "the %s '%s'" %(self.__class__.__name__, self.basename)
assert geoArr_mask.bands == 1, \
'Expected one single band as nodata mask for %s. Got %s bands.' % (self.basename, geoArr_mask.bands)
assert geoArr_mask.shape[:2] == self.shape[:2], 'The provided nodata mask must have the same number of ' \
'rows and columns as the %s itself.' %imName
assert geoArr_mask.gt == self.gt, \
'The geotransform of the given nodata mask for %s must match the geotransform of the %s itself. ' \
'Got %s.' %(imName, self.__class__.__name__, geoArr_mask.gt)
assert prj_equal(geoArr_mask.prj, self.prj), \
'The projection of the given nodata mask for the %s must match the projection of the %s itself.' \
%(imName, self.__class__.__name__)
self._mask_nodata = geoArr_mask
@property
......@@ -388,7 +410,7 @@ class GeoArray(object):
return self._footprint_poly
else:
assert self.mask_nodata is not None, 'A nodata mask is needed for calculating the footprint polygon. '
if np.std(self.mask_nodata)==0:
if np.std(self.mask_nodata[:])==0:
# do not run raster2polygon if whole image is filled with data
self._footprint_poly = self.box.mapPoly
else:
......@@ -526,6 +548,8 @@ class GeoArray(object):
if attr in self.__dir__(): #__dir__() includes also methods and properties
return self.__getattribute__(attr) #__getattribute__ avoids infinite loop
elif self.is_inmem and hasattr(np.array([]),attr):
return self[:].__getattribute__(attr)
else:
raise AttributeError("%s object has no attribute '%s'." %(self.__class__.__name__, attr))
......@@ -559,10 +583,10 @@ class GeoArray(object):
arr = self[:,:,fromBand] if self.ndim==3 else self[:]
if self.nodata is None:
self.mask_nodata = np.ones((self.rows, self.cols), np.uint8)
self.mask_nodata = np.ones((self.rows, self.cols), np.bool)
else:
self.mask_nodata = np.where(arr == self.nodata, 0, 1).astype(np.uint8) if arr.ndim == 2 else \
np.all(np.where(arr == self.nodata, 0, 1), axis=2).astype(np.uint8)
self.mask_nodata = np.where(arr == self.nodata, 0, 1).astype(np.bool) if arr.ndim == 2 else \
np.all(np.where(arr == self.nodata, 0, 1), axis=2).astype(np.bool)
def set_gdalDataset_meta(self):
......@@ -726,9 +750,9 @@ class GeoArray(object):
# set metadata
if not self.metadata.empty:
for bn, bidx in self.bandnames.items():
for bidx in range(self.bands):
band = ds.GetRasterBand(bidx+1)
band.SetMetadata(self.metadata[bn].to_dict())
band.SetMetadata(self.metadata[bidx].to_dict())
band = None
driver.CreateCopy(out_path, ds, options=creationOptions if creationOptions else [])
......@@ -1047,6 +1071,53 @@ class GeoArray(object):
return sub_arr, sub_gt, sub_prj
def reproject_to_new_grid(self, prototype=None, tgt_prj=None, tgt_xygrid=None, rspAlg='cubic'):
"""Reproject all array-like attributes to a given target grid.
:param prototype: <GeoArray> an instance of GeoArray to be used as pixel grid reference
:param tgt_prj: <str> WKT string of the projection
:param tgt_xygrid: <list> target XY grid, e.g. [[xmin,xmax], [ymax, ymin]] for the UL corner
:param rspAlg: <str, int> GDAL compatible resampling algorithm code
:return:
"""
assert (tgt_prj and tgt_xygrid) or prototype, "Provide either 'prototype' or 'tgt_prj' and 'tgt_xygrid'!"
tgt_prj = tgt_prj if tgt_prj else prototype.prj
tgt_xygrid = tgt_xygrid if tgt_xygrid is not None else prototype.xygrid_specs
assert tgt_xygrid[1][0]>tgt_xygrid[1][1]
# set target GSD
tgt_xgsd, tgt_ygsd = abs(tgt_xygrid[0][0]-tgt_xygrid[0][1]), abs(tgt_xygrid[1][0]-tgt_xygrid[1][1])
# set target bounds
tgt_bounds = reproject_shapelyGeometry(self.box.mapPoly, self.prj, tgt_prj).bounds
gt = (tgt_xygrid[0][0], tgt_xgsd, 0, max(tgt_xygrid[1]), 0, -tgt_ygsd)
xmin, ymin, xmax, ymax = snap_bounds_to_pixGrid(tgt_bounds, gt, roundAlg='on')
from ...geo.raster.reproject import warp_ndarray
self.arr, self.gt, self.prj = \
warp_ndarray(self[:], self.gt, self.prj, tgt_prj,
out_gsd = (tgt_xgsd, tgt_ygsd),
out_bounds = (xmin, ymin, xmax, ymax),
out_bounds_prj = tgt_prj,
rspAlg = rspAlg,
in_nodata = self.nodata,
CPUs = None,# if self.mp else 1, # TODO
progress = self.progress)
if hasattr(self, '_mask_nodata') and self._mask_nodata is not None:
self.mask_nodata.reproject_to_new_grid(prototype = prototype,
tgt_prj = tgt_prj,
tgt_xygrid = tgt_xygrid,
rspAlg = 'near')
if hasattr(self, '_mask_baddata') and self._mask_baddata is not None:
self.mask_baddata.reproject_to_new_grid(prototype = prototype,
tgt_prj = tgt_prj,
tgt_xygrid = tgt_xygrid,
rspAlg = 'near')
def read_pointData(self, mapXY_points, mapXY_points_prj=None, band=None):
"""Returns the array values for the given set of X/Y coordinates.
NOTE: If GeoArray has been instanced with a file path, the function will read the dataset into memory.
......@@ -1124,6 +1195,9 @@ class BadDataMask(GeoArray):
self.arr = self.arr.astype(np.bool)
# del self._mask_baddata, self.mask_baddata # TODO delete property (requires deleter)
@property
def arr(self):
return self._arr
......@@ -1135,13 +1209,46 @@ class BadDataMask(GeoArray):
assert len(pixelVals_in_mask) <= 2, 'Bad data mask must have only two pixel values (boolean) - 0 and 1 or ' \
'False and True! The given mask for %s contains the values %s.' % (
self.basename, pixelVals_in_mask)
assert pixelVals_in_mask in [[0, 1], [False, True]], 'Found unsupported pixel values in the given bad ' \
'data mask for %s: %s Only the values True, False, 0 ' \
'and 1 are supported. ' % (self.basename, pixelVals_in_mask)
assert pixelVals_in_mask in [[0, 1], [0],[1], [False, True], [False], [True]],\
'Found unsupported pixel values in the given bad data mask for %s: %s Only the values True, False, 0 ' \
'and 1 are supported. ' % (self.basename, pixelVals_in_mask)
self._arr = ndarray.astype(np.bool)
class NoDataMask(GeoArray):
def __init__(self, path_or_array, geotransform=None, projection=None, bandnames=None, nodata=None, progress=True,
q=False):
super(NoDataMask, self).__init__(path_or_array, geotransform=geotransform, projection=projection,
bandnames=bandnames, nodata=nodata, progress=progress, q=q)
self.arr = self.arr.astype(np.bool)
# del self._mask_nodata, self.mask_nodata # TODO delete property (requires deleter)
# TODO check that: "Automatically detected nodata value for NoDataMask 'IN_MEM': 1.0"
@property
def arr(self):
return self._arr
@arr.setter
def arr(self, ndarray):
assert isinstance(ndarray, np.ndarray), "'arr' can only be set to a numpy array!"
pixelVals_in_mask = sorted(list(np.unique(ndarray)))
assert len(pixelVals_in_mask) <= 2, 'Nodata mask must have only two pixel values (boolean) - 0 and 1 or ' \
'False and True! The given mask for %s contains the values %s.' % (
self.basename, pixelVals_in_mask)
assert pixelVals_in_mask in [[0, 1], [0], [1], [False, True], [False], [True]],\
'Found unsupported pixel values in the given Nodata mask for %s: %s Only the values True, False, 0 ' \
'and 1 are supported. ' % (self.basename, pixelVals_in_mask)
self._arr = ndarray.astype(np.bool)
def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0):
"""
NOTE: asserts that mapBounds have the same projection like the coordinates in arr_gt
......
......@@ -14,7 +14,7 @@ def find_nearest(array, value, roundAlg='auto', extrapolate=False, exclude_val=F
:param array:
:param value:
:param roundAlg:
:param roundAlg: <str> 'auto', 'on', 'off'
:param extrapolate: extrapolate the given array if the given value is outside the array
:param exclude_val:
"""
......
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