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

geo,map_info:

- geotransform2mapinfo(): bugfix

geo,projection:
- WKT2EPSG(): added type assertion
- modified assertion message

io.raster.GeoArray.GeoArray:
- arr.setter: better assertion message
- projection.setter: modified assertion (now more robust)
- mask_nodata(): bugfix for not setting projection properly
- mask_baddata(): bugfix for not setting projection properly
- added __setitem__() to overwrite pixel values with a given array
- __getattr__(): bugfix
- calc_mask_nodata(): bugfix
- set_gdalDataset_meta(): read projection is now extra validated
- from_path: added TODO
- _get_plottable_image(): bugfix
- to_disk(): now sets _arr directly

io.raster.GeoArray.BadDataMask:
- __init__(): bugfix

io.raster.GeoArray.NoDataMask:
- __init__(): bugfix; added TODO

- added class 'CloudMask'

- updated __version__
parent c235e9ca
...@@ -15,7 +15,7 @@ __all__=[#'compatibility', ...@@ -15,7 +15,7 @@ __all__=[#'compatibility',
'similarity', 'similarity',
'GeoArray'] 'GeoArray']
__version__ = '20170103_01' __version__ = '20170104_01'
__author__='Daniel Scheffler' __author__='Daniel Scheffler'
# Validate GDAL version # Validate GDAL version
......
...@@ -18,7 +18,7 @@ def geotransform2mapinfo(gt,prj): ...@@ -18,7 +18,7 @@ def geotransform2mapinfo(gt,prj):
:rtype: list :rtype: list
""" """
if gt is None: if gt in [None, [0, 1, 0, 0, 0, -1]]:
return None return None
if gt[2]!=0 or gt[4]!=0: # TODO if gt[2]!=0 or gt[4]!=0: # TODO
......
...@@ -96,11 +96,14 @@ def WKT2EPSG(wkt, epsg=os.environ['GDAL_DATA'].replace('/gdal', '/proj/epsg')): ...@@ -96,11 +96,14 @@ def WKT2EPSG(wkt, epsg=os.environ['GDAL_DATA'].replace('/gdal', '/proj/epsg')):
""" """
# FIXME this function returns None if datum=NAD27 but works with datum=WGS84, e.g.: # FIXME this function returns None if datum=NAD27 but works with datum=WGS84, e.g.:
# FIXME {PROJCS["UTM_Zone_33N",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]} # FIXME {PROJCS["UTM_Zone_33N",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]}
if not isinstance(wkt,str):
raise TypeError("'wkt' must be a string. Received %s." %type(wkt))
code = None #default code = None #default
p_in = osr.SpatialReference() p_in = osr.SpatialReference()
s = p_in.ImportFromWkt(wkt) s = p_in.ImportFromWkt(wkt)
if s == 5: # invalid WKT if s == 5: # invalid WKT
raise Exception('Invalid WKT.') raise Exception('Received an invalid WKT string: %s' %wkt)
if p_in.IsLocal(): if p_in.IsLocal():
raise Exception('The given WKT is a local coordinate system.') raise Exception('The given WKT is a local coordinate system.')
cstype = 'GEOGCS' if p_in.IsGeographic() else 'PROJCS' cstype = 'GEOGCS' if p_in.IsGeographic() else 'PROJCS'
......
...@@ -132,7 +132,7 @@ class GeoArray(object): ...@@ -132,7 +132,7 @@ class GeoArray(object):
@arr.setter @arr.setter
def arr(self, ndarray): def arr(self, ndarray):
assert isinstance(ndarray, np.ndarray), "'arr' can only be set to a numpy array!" assert isinstance(ndarray, np.ndarray), "'arr' can only be set to a numpy array! Got %s." %type(ndarray)
self._arr = ndarray self._arr = ndarray
...@@ -274,8 +274,9 @@ class GeoArray(object): ...@@ -274,8 +274,9 @@ class GeoArray(object):
@projection.setter @projection.setter
def projection(self, prj): def projection(self, prj):
if self.filePath: if self.filePath:
assert self.projection in [None,prj], "Cannot set %s.projection to the given value because it does not " \ assert self.projection is None or prj_equal(self.projection, prj),\
"match the projection from the file on disk." %self.__class__.__name__ "Cannot set %s.projection to the given value because it does not match the projection from the file " \
"on disk." %self.__class__.__name__
else: else:
self._projection = prj self._projection = prj
...@@ -349,7 +350,7 @@ class GeoArray(object): ...@@ -349,7 +350,7 @@ class GeoArray(object):
if mask is not None: if mask is not None:
geoArr_mask = NoDataMask(mask, progress=self.progress, q=self.q) geoArr_mask = NoDataMask(mask, progress=self.progress, q=self.q)
geoArr_mask.gt = geoArr_mask.gt if geoArr_mask.gt not in [None, [0, 1, 0, 0, 0, -1]] else self.gt 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 geoArr_mask.prj = geoArr_mask.prj if geoArr_mask.prj else self.prj
imName = "the %s '%s'" %(self.__class__.__name__, self.basename) imName = "the %s '%s'" %(self.__class__.__name__, self.basename)
assert geoArr_mask.bands == 1, \ assert geoArr_mask.bands == 1, \
...@@ -359,7 +360,7 @@ class GeoArray(object): ...@@ -359,7 +360,7 @@ class GeoArray(object):
assert geoArr_mask.gt == self.gt, \ assert geoArr_mask.gt == self.gt, \
'The geotransform of the given nodata mask for %s must match the geotransform of the %s itself. ' \ '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) 'Got %s.' %(imName, self.__class__.__name__, geoArr_mask.gt)
assert prj_equal(geoArr_mask.prj, self.prj), \ assert not geoArr_mask.prj or 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.' \ 'The projection of the given nodata mask for the %s must match the projection of the %s itself.' \
%(imName, self.__class__.__name__) %(imName, self.__class__.__name__)
...@@ -384,7 +385,7 @@ class GeoArray(object): ...@@ -384,7 +385,7 @@ class GeoArray(object):
if mask is not None: if mask is not None:
geoArr_mask = BadDataMask(mask, progress=self.progress, q=self.q) geoArr_mask = BadDataMask(mask, progress=self.progress, q=self.q)
geoArr_mask.gt = geoArr_mask.gt if geoArr_mask.gt not in [None, [0, 1, 0, 0, 0, -1]] else self.gt 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 geoArr_mask.prj = geoArr_mask.prj if geoArr_mask.prj else self.prj
imName = "the %s '%s'" %(self.__class__.__name__, self.basename) imName = "the %s '%s'" %(self.__class__.__name__, self.basename)
assert geoArr_mask.bands == 1, \ assert geoArr_mask.bands == 1, \
...@@ -541,6 +542,20 @@ class GeoArray(object): ...@@ -541,6 +542,20 @@ class GeoArray(object):
return self._arr_cache return self._arr_cache
def __setitem__(self, idx, array2set):
"""Overwrites the pixel values of GeoArray.arr with the given array.
:param idx: <int, list, slice> the index position to overwrite
:param array2set: <np.ndarray> array to be set. Must be compatible to the given index position.
:return:
"""
if self.is_inmem:
self.arr[idx] = array2set
else:
raise NotImplementedError('Item assignment for %s instances that are not in memory is not yet supported.'
%self.__class__.__name__)
def __getattr__(self, attr): def __getattr__(self, attr):
# check if the requested attribute can not be present because GeoArray has been instanced with an array # check if the requested attribute can not be present because GeoArray has been instanced with an array
if attr not in self.__dir__() and not self.is_inmem and attr in ['shape','dtype','geotransform', 'projection']: if attr not in self.__dir__() and not self.is_inmem and attr in ['shape','dtype','geotransform', 'projection']:
...@@ -548,7 +563,7 @@ class GeoArray(object): ...@@ -548,7 +563,7 @@ class GeoArray(object):
if attr in self.__dir__(): #__dir__() includes also methods and properties if attr in self.__dir__(): #__dir__() includes also methods and properties
return self.__getattribute__(attr) #__getattribute__ avoids infinite loop return self.__getattribute__(attr) #__getattribute__ avoids infinite loop
elif self.is_inmem and hasattr(np.array([]),attr): elif hasattr(np.array([]),attr):
return self[:].__getattribute__(attr) return self[:].__getattribute__(attr)
else: else:
raise AttributeError("%s object has no attribute '%s'." %(self.__class__.__name__, attr)) raise AttributeError("%s object has no attribute '%s'." %(self.__class__.__name__, attr))
...@@ -580,7 +595,7 @@ class GeoArray(object): ...@@ -580,7 +595,7 @@ class GeoArray(object):
""" """
if self._mask_nodata is None or overwrite: if self._mask_nodata is None or overwrite:
assert self.ndim in [2, 3], "Only 2D or 3D arrays are supported. Got a %sD array." % self.ndim assert self.ndim in [2, 3], "Only 2D or 3D arrays are supported. Got a %sD array." % self.ndim
arr = self[:,:,fromBand] if self.ndim==3 else self[:] arr = self[:,:,fromBand] if self.ndim==3 and fromBand is not None else self[:]
if self.nodata is None: if self.nodata is None:
self.mask_nodata = np.ones((self.rows, self.cols), np.bool) self.mask_nodata = np.ones((self.rows, self.cols), np.bool)
...@@ -602,7 +617,7 @@ class GeoArray(object): ...@@ -602,7 +617,7 @@ class GeoArray(object):
self._shape = tuple([ds.RasterYSize, ds.RasterXSize] + ([ds.RasterCount] if ds.RasterCount>1 else [])) self._shape = tuple([ds.RasterYSize, ds.RasterXSize] + ([ds.RasterCount] if ds.RasterCount>1 else []))
self._dtype = gdal_array.GDALTypeCodeToNumericTypeCode(ds.GetRasterBand(1).DataType) self._dtype = gdal_array.GDALTypeCodeToNumericTypeCode(ds.GetRasterBand(1).DataType)
self._geotransform = ds.GetGeoTransform() self._geotransform = ds.GetGeoTransform()
self._projection = ds.GetProjection() self._projection = EPSG2WKT(WKT2EPSG(ds.GetProjection())) # temp conversion to EPSG needed because GDAL seems to modify WKT string when writing file to disk (e.g. using gdal_merge) -> conversion to EPSG and back undos that
if not 'nodata' in self._initParams or self._initParams['nodata'] is None: if not 'nodata' in self._initParams or self._initParams['nodata'] is None:
band = ds.GetRasterBand(1) band = ds.GetRasterBand(1)
self._nodata = band.GetNoDataValue() # FIXME this does not support different nodata values within the same file self._nodata = band.GetNoDataValue() # FIXME this does not support different nodata values within the same file
...@@ -723,7 +738,8 @@ class GeoArray(object): ...@@ -723,7 +738,8 @@ class GeoArray(object):
if out_arr.shape==self.shape: if out_arr.shape==self.shape:
self.arr = out_arr self.arr = out_arr
return out_arr return out_arr # TODO implement check of returned datatype (e.g. NoDataMask should always return np.bool
# TODO -> would be np.int8 if an int8 file is read from disk
def save(self, out_path, fmt='ENVI', creationOptions=None): def save(self, out_path, fmt='ENVI', creationOptions=None):
...@@ -793,7 +809,7 @@ class GeoArray(object): ...@@ -793,7 +809,7 @@ class GeoArray(object):
cS, cE = xlim if isinstance(xlim, (tuple, list)) else (0, self.columns - 1) 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) rS, rE = ylim if isinstance(ylim, (tuple, list)) else (0, self.rows - 1)
image2plot = self[rS:rE, cS:cE, band] image2plot = self[rS:rE, cS:cE, band] if band is not None else self[rS:rE, cS:cE]
gt, prj = self.geotransform, self.projection gt, prj = self.geotransform, self.projection
...@@ -871,8 +887,9 @@ class GeoArray(object): ...@@ -871,8 +887,9 @@ class GeoArray(object):
rS, rE = ylim if isinstance(ylim, (tuple, list)) else (0, self.rows - 1) rS, rE = ylim if isinstance(ylim, (tuple, list)) else (0, self.rows - 1)
image2plot = np.array(rescale_intensity(image2plot, in_range=(vmin, vmax))) image2plot = np.array(rescale_intensity(image2plot, in_range=(vmin, vmax)))
get_hv_image = lambda b: hv.Image(image2plot[:,:,b], bounds=(cS, rS, cE, rE))(style={'cmap': 'gray'}, # FIXME ylabels have the wrong order get_hv_image = lambda b: hv.Image(image2plot[:,:,b] if b is not None else image2plot,
plot={'fig_inches':4 if figsize is None else figsize, 'show_grid':True}) bounds=(cS, rS, cE, rE))(style={'cmap': 'gray'}, # FIXME ylabels have the wrong order
plot={'fig_inches':4 if figsize is None else figsize, 'show_grid':True})
#hvIm = hv.Image(image2plot)(style={'cmap': 'gray'}, figure_inches=figsize) #hvIm = hv.Image(image2plot)(style={'cmap': 'gray'}, figure_inches=figsize)
hmap = hv.HoloMap([(band, get_hv_image(band)) for band in range(image2plot.shape[2])], kdims=['band']) hmap = hv.HoloMap([(band, get_hv_image(band)) for band in range(image2plot.shape[2])], kdims=['band'])
...@@ -1167,12 +1184,12 @@ class GeoArray(object): ...@@ -1167,12 +1184,12 @@ class GeoArray(object):
return self return self
def _to_disk(self): def to_disk(self):
"""Sets self.arr back to None if GeoArray has been instanced with a file path """Sets self.arr back to None if GeoArray has been instanced with a file path
and the whole datadaset has been read.""" and the whole dataset has been read."""
if self.filePath and os.path.isfile(self.filePath): if self.filePath and os.path.isfile(self.filePath):
self.arr = None self._arr = None
else: else:
warnings.warn('GeoArray object cannot be turned into disk mode because this asserts that GeoArray.filePath ' warnings.warn('GeoArray object cannot be turned into disk mode because this asserts that GeoArray.filePath '
'contains a valid file path. Got %s.' %self.filePath) 'contains a valid file path. Got %s.' %self.filePath)
...@@ -1206,7 +1223,8 @@ class BadDataMask(GeoArray): ...@@ -1206,7 +1223,8 @@ class BadDataMask(GeoArray):
super(BadDataMask, self).__init__(path_or_array, geotransform=geotransform, projection=projection, super(BadDataMask, self).__init__(path_or_array, geotransform=geotransform, projection=projection,
bandnames=bandnames, nodata=nodata, progress=progress, q=q) bandnames=bandnames, nodata=nodata, progress=progress, q=q)
self.arr = self.arr.astype(np.bool) if self.is_inmem:
self.arr = self.arr.astype(np.bool)
# del self._mask_baddata, self.mask_baddata # TODO delete property (requires deleter) # del self._mask_baddata, self.mask_baddata # TODO delete property (requires deleter)
...@@ -1237,11 +1255,15 @@ class NoDataMask(GeoArray): ...@@ -1237,11 +1255,15 @@ class NoDataMask(GeoArray):
super(NoDataMask, self).__init__(path_or_array, geotransform=geotransform, projection=projection, super(NoDataMask, self).__init__(path_or_array, geotransform=geotransform, projection=projection,
bandnames=bandnames, nodata=nodata, progress=progress, q=q) bandnames=bandnames, nodata=nodata, progress=progress, q=q)
self.arr = self.arr.astype(np.bool) if self.is_inmem:
self.arr = self.arr.astype(np.bool)
# del self._mask_nodata, self.mask_nodata # TODO delete property (requires deleter) # 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" # TODO check that: "Automatically detected nodata value for NoDataMask 'IN_MEM': 1.0"
# TODO disk-mode: init must check the numbers of bands, and ideally also the pixel values in mask
# TODO nodata=False?
@property @property
def arr(self): def arr(self):
...@@ -1261,6 +1283,20 @@ class NoDataMask(GeoArray): ...@@ -1261,6 +1283,20 @@ class NoDataMask(GeoArray):
class CloudMask(GeoArray):
def __init__(self, path_or_array, geotransform=None, projection=None, bandnames=None, nodata=None, progress=True,
q=False):
# TODO implement class definitions and specific metadata
super(CloudMask, self).__init__(path_or_array, geotransform=geotransform, projection=projection,
bandnames=bandnames, nodata=nodata, progress=progress, q=q)
# del self._mask_nodata, self.mask_nodata # TODO delete property (requires deleter)
# TODO check that: "Automatically detected nodata value for CloudMask 'IN_MEM': 1.0"
def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0): def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0):
""" """
......
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