Commit 42ae52b1 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

simplified GeoArray.__init__(); some bugfixes and further developments

geo.coord_grid:
- find_nearest_grid_coord(): fixed inconsistency regarding key 'SE' vs. 'SO'

io.raster.GeoArray.GeoArray:
- __init__():
    - added docstring
    - init parameters are now respected in case GeoArray is instanced from another GeoArray
    - moved some parameter settings to separate properties
- attributes 'arr' and 'bandnames' are now properties
- revised propery 'shape'
- show(): added warning
- added deepcopy_array()
parent 44f166b7
......@@ -44,8 +44,8 @@ def is_point_on_grid(pointXY,xgrid,ygrid):
def find_nearest_grid_coord(valXY,gt,rows,cols,direction='NW',extrapolate=True):
UL, LL, LR, UR = get_corner_coordinates(gt=gt, rows=rows, cols=cols) # (x,y) tuples
round_x = {'NW': 'off', 'NO': 'on', 'SW': 'off', 'SO': 'on'}[direction]
round_y = {'NW': 'on', 'NO': 'on', 'SW': 'off', 'SO': 'off'}[direction]
round_x = {'NW': 'off', 'NO': 'on', 'SW': 'off', 'SE': 'on' }[direction]
round_y = {'NW': 'on', 'NO': 'on', 'SW': 'off', 'SE': 'off'}[direction]
tgt_xgrid = np.arange(UL[0], UR[0] + gt[1], gt[1])
tgt_ygrid = np.arange(LL[1], UL[1] + abs(gt[5]), abs(gt[5]))
tgt_x = find_nearest(tgt_xgrid, valXY[0], roundAlg=round_x, extrapolate=extrapolate)
......@@ -53,7 +53,7 @@ def find_nearest_grid_coord(valXY,gt,rows,cols,direction='NW',extrapolate=True):
return tgt_x,tgt_y
def move_shapelyPoly_to_image_grid(shapelyPoly, gt, rows, cols, moving_dir='NW'):
def move_shapelyPoly_to_image_grid(shapelyPoly, gt, rows, cols, moving_dir='NW', extrapolate=True):
polyULxy = (min(shapelyPoly.exterior.coords.xy[0]), max(shapelyPoly.exterior.coords.xy[1]))
polyLRxy = (max(shapelyPoly.exterior.coords.xy[0]), min(shapelyPoly.exterior.coords.xy[1]))
UL, LL, LR, UR = get_corner_coordinates(gt=gt, rows=rows, cols=cols) # (x,y) tuples
......@@ -61,8 +61,8 @@ def move_shapelyPoly_to_image_grid(shapelyPoly, gt, rows, cols, moving_dir='NW')
round_y = {'NW': 'on', 'NO': 'on', 'SW': 'off', 'SE': 'off'}[moving_dir]
tgt_xgrid = np.arange(UL[0], UR[0] + gt[1], gt[1])
tgt_ygrid = np.arange(LL[1], UL[1] + abs(gt[5]), abs(gt[5]))
tgt_xmin = find_nearest(tgt_xgrid, polyULxy[0], roundAlg=round_x)
tgt_xmax = find_nearest(tgt_xgrid, polyLRxy[0], roundAlg=round_x)
tgt_ymin = find_nearest(tgt_ygrid, polyLRxy[1], roundAlg=round_y)
tgt_ymax = find_nearest(tgt_ygrid, polyULxy[1], roundAlg=round_y)
tgt_xmin = find_nearest(tgt_xgrid, polyULxy[0], roundAlg=round_x, extrapolate=True)
tgt_xmax = find_nearest(tgt_xgrid, polyLRxy[0], roundAlg=round_x, extrapolate=True)
tgt_ymin = find_nearest(tgt_ygrid, polyLRxy[1], roundAlg=round_y, extrapolate=True)
tgt_ymax = find_nearest(tgt_ygrid, polyULxy[1], roundAlg=round_y, extrapolate=True)
return box(tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax)
\ No newline at end of file
......@@ -49,8 +49,12 @@ def _alias_property(key):
class GeoArray(object):
def __init__(self, path_or_array, geotransform=None, projection=None, bandnames=None, nodata=None, progress=True,
q=False):
# type: (Any, tuple, str, list) -> GeoArray
"""
# type: (Any, tuple, str, list, float, bool, bool) -> GeoArray
"""This class creates a fast Python inteface to geodata - either on disk or in memory. It can be instanced with
a file path or with a numpy array and the corresponding geoinformation. Instances can always be indexed like
normal numpy arrays, no matter if GeoArray has been instanced from file or from an in-memory array. GeoArray
provides a wide range of geo-related attributes belonging to the dataset as well as some functions for quickly
visualizing the data as a map, a simple image or an interactive image.
:param path_or_array: a numpy.ndarray or a valid file path
:param geotransform: GDAL geotransform of the given array or file on disk
......@@ -66,17 +70,32 @@ class GeoArray(object):
# FIXME implement compatibility to GDAL VRTs
if type(path_or_array) not in [str, np.ndarray, type(self)]:
if type(path_or_array) not in [str, np.ndarray, type(GeoArray)]:
raise ValueError("%s parameter 'arg' takes only string "
"or np.ndarray types. Got %s." %(self.__class__.__name__,type(path_or_array)))
if isinstance(path_or_array, str):
assert ' ' not in path_or_array, "The given path contains whitespaces. This is not supported by GDAL."
if not os.path.exists(path_or_array):
raise FileNotFoundError(path_or_array) if PY3 else FileNotFoundError_comp(path_or_array)
if type(path_or_array)==type(self):
self.__dict__= path_or_array.__dict__
self._initParams = dict([x for x in locals().items() if x[0] != "self"])
self.geotransform = geotransform if geotransform else self.geotransform
self.projection = projection if projection else self.projection
self.bandnames = bandnames if bandnames else self.bandnames
self._nodata = nodata if nodata is not None else self._nodata
self.progress = progress if progress else self.progress
self.q = q if q is not None else self.q
else:
self._initParams = dict([x for x in locals().items() if x[0] != "self"])
self.arg = path_or_array
self.arr = self.arg if isinstance(self.arg, np.ndarray) else None
self.filePath = self.arg if isinstance(self.arg, str) and self.arg else None
self._arr = path_or_array if isinstance(path_or_array, np.ndarray) else None
self.filePath = path_or_array if isinstance(path_or_array, str) and path_or_array else None
self.basename = os.path.splitext(os.path.basename(self.filePath))[0] if not self.is_inmem else 'IN_MEM'
self.progress = progress
self.q = q
......@@ -90,41 +109,48 @@ class GeoArray(object):
self._footprint_poly = None
self._gdalDataset_meta_already_set = False
self._metadata = None
if isinstance(self.arg, str):
assert ' ' not in self.arg, "The given path contains whitespaces. This is not supported by GDAL."
if not os.path.exists(self.filePath):
raise FileNotFoundError(self.filePath) if PY3 else FileNotFoundError_comp(self.filePath)
ds = gdal.Open(self.filePath)
if not ds:
raise Exception('Error reading file: ' + gdal.GetLastErrorMsg())
bands = ds.RasterCount
ds = None
else:
self._shape = self.arr.shape
self._dtype = self.arr.dtype
bands = self.arr.shape[2] if len(self.arr.shape) == 3 else 1
self._bandnames = None
if bandnames:
assert len(bandnames) == bands, \
'Number of given bandnames does not match number of bands in array.'
assert len(list(set([type(b) for b in bandnames]))) == 1 and type(bandnames[0] == 'str'), \
"'bandnames must be a set of strings. Got other datetypes in there.'"
self.bandnames = {band: i for i, band in enumerate(bandnames)} # syntax supported since Python 2.7
assert len(self.bandnames) == bands, 'Bands must not have the same name.'
else:
self.bandnames = {'B%s' % band: i for i, band in enumerate(range(1, bands + 1))}
self.bandnames = bandnames # use property in order to validate given value
if geotransform:
self.geotransform = geotransform # use property in order to validate given value
if projection:
self.projection = projection # use property in order to validate given value
self.projection = projection # use property in order to validate given value
if self.filePath:
self.set_gdalDataset_meta()
@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!"
self._arr = ndarray
@property
def bandnames(self):
if self._bandnames:
return self._bandnames
else:
self._bandnames = {'B%s' % band: i for i, band in enumerate(range(1, self.bands + 1))}
return self._bandnames
@bandnames.setter
def bandnames(self, list_bandnames):
assert len(list_bandnames) == self.bands, \
'Number of given bandnames does not match number of bands in array.'
assert len(list(set([type(b) for b in list_bandnames]))) == 1 and type(list_bandnames[0] == 'str'), \
"'bandnames must be a set of strings. Got other datetypes in there.'"
bN_dict = {band: i for i, band in enumerate(list_bandnames)} # syntax supported since Python 2.7
assert len(bN_dict) == self.bands, 'Bands must not have the same name.'
self._bandnames = bN_dict
@property
def is_inmem(self):
......@@ -135,11 +161,14 @@ class GeoArray(object):
@property
def shape(self):
"""Get the array shape of the associated image array."""
if self._shape:
return self._shape
if self.is_inmem:
return self.arr.shape
else:
self.set_gdalDataset_meta()
return self._shape
if self._shape:
return self._shape
else:
self.set_gdalDataset_meta()
return self._shape
@property
......@@ -312,6 +341,8 @@ class GeoArray(object):
@property
def footprint_poly(self):
# FIXME should return polygon in image coordiates if no projection is available
"""Get the footprint polygon of the associated image array (returns an instance of shapely.geometry.Polygon."""
if self._footprint_poly is not None:
return self._footprint_poly
......@@ -782,6 +813,10 @@ class GeoArray(object):
return hmap
else:
if interactive:
warnings.warn('Currently there is no interactive mode for single-band arrays. '
'Switching to standard matplotlib figure..') # TODO implement zoomable fig
# show image
plt.figure(figsize=figsize)
rows, cols = image2plot.shape[:2]
......@@ -1020,6 +1055,13 @@ class GeoArray(object):
return self
def deepcopy_array(self):
if self.is_inmem:
temp = np.empty_like(self.arr)
temp[:] = self.arr
self.arr = temp # deep copy: converts view to its own array in order to avoid wrong output
def cache_array_subset(self, subarray):
"""Sets the array cache of the GeoArray instance to the given array in order to speed up calculations
afterwards."""
......
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