added warning if automatic nodata value detection returns unreliable value;...

added warning if automatic nodata value detection returns unreliable value; bugfix for overwriting user defined attributes of GeoArray

io.raster.GeoArray:
- GeoArray:
    - added quiet mode
    - added attribute _gdalDataset_meta_already_set
    - added print output if nodata value had to be detected automatically
    - set_gdalDataset_meta(): added docstring; method is now only called from __init__

numeric.array:
- revised find_noDataVal()

- added new package similarity

ptds/__init__.py

-#import py_tools_ds.ptds.io
-#import py_tools_ds.ptds.geo
-#import py_tools_ds.ptds.numeric
-#import py_tools_ds.ptds.processing
+#from . import compatibility
 from . import io
 from . import geo
 from . import numeric
 from . import processing
+from . import similarity
 
 from .io.raster.GeoArray import GeoArray
-#from py_tools_ds.ptds.io.raster.GeoArray import GeoArray
 
-__all__=['io',
+__all__=[#'compatibility',
+         'io',
          'geo',
          'numeric',
          'processing',
+         'similarity',
          'GeoArray']
 
 __version__ = '20161102_01'

ptds/io/raster/GeoArray.py

@@ -41,7 +41,7 @@ def _alias_property(key):
 
 
 class GeoArray(object):
-    def __init__(self, path_or_array, geotransform=None, projection=None, bandnames=None, nodata=None):
+    def __init__(self, path_or_array, geotransform=None, projection=None, bandnames=None, nodata=None, q=False):
         # type: (Any, tuple, str, list) -> GeoArray
         """
 
@@ -51,7 +51,8 @@ class GeoArray(object):
                                 (only needed if GeoArray is instanced with an array)
         :param bandnames:       names of the bands within the input array, e.g. ['mask_1bit', 'mask_clouds'],
                                 (default: ['B1', 'B2', 'B3', ...])
-        :param nodata           nodata value
+        :param nodata:          nodata value
+        :param q:               quiet mode (default: False)
         """
 
         # FIXME implement compatibility to GDAL VRTs
@@ -64,6 +65,7 @@ class GeoArray(object):
         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.basename        = os.path.splitext(os.path.basename(self.filePath))[0] if not self.is_inmem else 'IN_MEM'
+        self.q               = q
         self._arr_cache      = None
         self._geotransform   = None
         self._projection     = None
@@ -72,6 +74,7 @@ class GeoArray(object):
         self._nodata         = nodata
         self._mask_nodata    = None
         self._footprint_poly = None
+        self._gdalDataset_meta_already_set = False
 
         if isinstance(self.arg, str):
             assert ' ' not in self.arg, "The given path contains whitespaces. This is not supported by GDAL."
@@ -92,6 +95,9 @@ class GeoArray(object):
             self._dtype = self.arr.dtype
             bands = self.arr.shape[2] if len(self.arr.shape) == 3 else 1
 
+        if self.filePath:
+            self.set_gdalDataset_meta()
+
         if bandnames:
             assert len(bandnames) == bands, \
                 'Number of given bandnames does not match number of bands in array.'
@@ -172,11 +178,7 @@ class GeoArray(object):
         assert isinstance(gt,(list,tuple)) and len(gt)==6, 'geotransform must be a list with 6 numbers. Got %s.' %gt
         for i in gt: assert isinstance(i,(int,float)),     "geotransform must contain only numbers. Got '%s'." %i
 
-        if self.filePath:
-            assert self.geotransform == gt, "Cannot set %s.geotransform to the given value because it does not " \
-                                            "match the geotransform from the file on disk." %self.__class__.__name__
-        else:
-            self._geotransform = gt
+        self._geotransform = gt
 
 
     gt = _alias_property('geotransform')
@@ -244,6 +246,10 @@ class GeoArray(object):
                 if self._nodata == 'ambiguous':
                     warnings.warn('Nodata value could not be clearly identified. It has been set to None.')
                     self._nodata = None
+                else:
+                    if not self.q:
+                        print("Automatically detected nodata value for %s '%s': %s"
+                              %(self.__class__.__name__, self.basename, self._nodata))
             return self._nodata
 
 
@@ -368,16 +374,24 @@ class GeoArray(object):
 
 
     def set_gdalDataset_meta(self):
-        assert self.filePath
-        ds = gdal.Open(self.filePath)
-        # set private class variables (in order to avoid recursion error)
-        self._shape        = tuple([ds.RasterYSize, ds.RasterXSize] + ([ds.RasterCount] if ds.RasterCount>1 else []))
-        self._dtype        = gdal_array.GDALTypeCodeToNumericTypeCode(ds.GetRasterBand(1).DataType)
-        self._geotransform = ds.GetGeoTransform()
-        self._projection   = ds.GetProjection()
-        band               = ds.GetRasterBand(1)
-        self._nodata       = band.GetNoDataValue() # FIXME this does not support different nodata values within the same file
-        ds = band = None
+        """Retrieves GDAL metadata from file. This function is only executed once to avoid overwriting of user defined
+         attributes, that are defined after object instanciation.
+
+        :return:
+        """
+        if not self._gdalDataset_meta_already_set:
+            assert self.filePath
+            ds = gdal.Open(self.filePath)
+            # set private class variables (in order to avoid recursion error)
+            self._shape        = tuple([ds.RasterYSize, ds.RasterXSize] + ([ds.RasterCount] if ds.RasterCount>1 else []))
+            self._dtype        = gdal_array.GDALTypeCodeToNumericTypeCode(ds.GetRasterBand(1).DataType)
+            self._geotransform = ds.GetGeoTransform()
+            self._projection   = ds.GetProjection()
+            band               = ds.GetRasterBand(1)
+            self._nodata       = band.GetNoDataValue() # FIXME this does not support different nodata values within the same file
+            ds = band = None
+
+        self._gdalDataset_meta_already_set = True
 
 
     def from_path(self, path, getitem_params=None):
@@ -935,8 +949,8 @@ def get_array_at_mapPosOLD(arr, arr_gt, arr_prj, mapBounds, mapBounds_prj, band2
     return out_arr, out_gt, out_prj
 
 
-def get_array_at_mapPos(arr, arr_gt, arr_prj, out_prj, mapBounds, mapBounds_prj=None, out_gsd=None, band2get=None, fillVal=0,
-                        rspAlg='near'):
+def get_array_at_mapPos(arr, arr_gt, arr_prj, out_prj, mapBounds, mapBounds_prj=None, out_gsd=None, band2get=None,
+                        fillVal=0, rspAlg='near'):
     """
 
     :param arr:

ptds/numeric/array.py

@@ -2,6 +2,8 @@
 __author__ = "Daniel Scheffler"
 
 import numpy as np
+import warnings
+
 
 
 def get_outFillZeroSaturated(dtype):
@@ -28,12 +30,27 @@ def find_noDataVal(pathIm_or_GeoArray,bandIdx=0,sz=3):
     from .. import GeoArray
     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]
+
+    wins         = [geoArr[0:sz,0:sz,bandIdx], geoArr[0:sz,-sz:,bandIdx],
+                    geoArr[-sz:,-sz:,bandIdx], geoArr[-sz:,0:sz,bandIdx]] # UL, UR, LR, LL
+    means_stds   = [get_mean_std(win) for win in wins]
+
+    possVals     = [i['mean'] for i in means_stds 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 'ambiguous'
+
+    if possVals:
+        if np.std(possVals)!=0:
+            # different possible nodata values have been found in the image corner
+            return 'ambiguous'
+        else:
+            if len(possVals)<=2:
+                # each window in each corner
+                warnings.warn("\nAutomatic nodata value detection returned the value %s for GeoArray '%s' but this "
+                              "seems to be unreliable (occurs in only %s). To avoid automatic detection, just pass "
+                              "the correct nodata value." %(possVals[0], geoArr.basename,
+                                                            ('2 image corners' if len(possVals)==2 else '1 image corner')))
+            return possVals[0]
+    else:
+        return None
 
 

ptds/similarity/__init__.py

+# -*- coding: utf-8 -*-
+__author__ = "Daniel Scheffler"
+
+from . import raster
+
+__all__=['raster']
\ No newline at end of file

ptds/similarity/raster.py

+# -*- coding: utf-8 -*-
+__author__ = "Daniel Scheffler"
+
+from skimage.measure import compare_ssim as ssim
+
+
+def calc_ssim(match, other, dynamic_range=None, win_size=None):
+    """Calculates Mean Structural Similarity Index between two images.
+
+    :param match:
+    :param other:
+    :param dynamic_range:
+    :param win_size:
+    :return:
+    """
+    return ssim(match, other, dynamic_range=dynamic_range, win_size=win_size)
\ No newline at end of file