Merge branch 'dev'

ymraster/ymraster.py

@@ -50,6 +50,35 @@ import os
 import shutil
 from tempfile import gettempdir
 
+def add_suffix(path_file, suffix):
+    """
+    Parameters
+    ----------
+    path_file : str
+        A file path
+    suffix : str
+        Suffix to add at the end of a filename
+
+    Returns
+    -------
+    path_file_suffix : str
+        Path file with suffix add
+
+    Example
+    -------
+        >>> suffix = 'suffix'
+        >>> path_file = 'tata/toto.tif'
+        >>> suffix = 'suffix'
+        >>> path_file_suffix = add_suffix(path_file, suffix)
+        >>> print path_file
+        tata/toto.tif
+        >>> print path_file_suffix
+        tata/totosuffix.tif
+    """
+    begin, ext = os.path.splitext(path_file)
+    path_file_suffix = begin + suffix + ext
+
+    return path_file_suffix
 
 def _dt2float(dt):
     """Returns a float corresponding to the given datetime object.
@@ -259,8 +288,10 @@ def concatenate_rasters(rasters, **kw):
     srs, extent, otb_dtype = (raster0.srs,
                               raster0.gdal_extent,
                               raster0.dtype.otb_dtype)
-    if not srs :
-        print "Warning : Image has no Coordinate Reference System: '{:f}'".format(raster0)
+    if not srs:
+        warning_msg = ("Warning : Image has no Coordinate Reference System:" +
+                       "'{:f}'".format(raster0))
+        print warning_msg
 
     same_type = True
     for raster in rasters:
@@ -272,23 +303,25 @@ def concatenate_rasters(rasters, **kw):
             "Images have not the same extent: "
         "'{:f}' and '{:f}'".format(raster0, raster)
         if raster.dtype.otb_dtype != otb_dtype:
-            print "Warnig: Images have not same data type, float data type is assigned by default"
+            warning_msg = ('Warning: Images have not same data type, float' +
+                           ' data type is assigned by default')
+            print warning_msg
             same_type = False
 
     # Out file
-    out_filename = kw['out_filename'] \
-        if kw.get('out_filename') \
-        else os.path.join(gettempdir(), 'concat.tif')
+    out_filename = kw.get('out_filename',
+                          os.path.join(gettempdir(), 'concat.tif'))
 
-    #remove bands if indicated
+    # remove bands if indicated
     filenames = []
     bands_to_keep = kw.get('bands', None)
     if bands_to_keep:
-        for i,(band, rst) in enumerate(zip(bands_to_keep,rasters)):
+        for i, (band, rst) in enumerate(zip(bands_to_keep, rasters)):
             if band:
-                out = os.path.join(gettempdir(),'rst_reordered_{}.tif'.format(i+1))
+                out = os.path.join(gettempdir(),
+                                   'rst_reordered_{}.tif'.format(i+1))
                 filenames.append(out)
-                rst.extract_band(*band, out_filename = out)
+                rst.extract_band(*band, out_filename=out)
 
     # Perform the concatenation
     if not filenames:
@@ -302,7 +335,7 @@ def concatenate_rasters(rasters, **kw):
         ConcatenateImages.SetParameterOutputImagePixelType("out", otb_dtype)
     ConcatenateImages.ExecuteAndWriteOutput()
 
-    #Remove temp file if needed
+    # Remove temp file if needed
     if bands_to_keep:
         for f in filenames:
             os.remove(f)
@@ -604,6 +637,32 @@ class Raster(Sized):
         # Close file
         ds = None
 
+    def band_statistics(self, idx):
+        """
+        Compute band statistic, without taking into account nodata values
+
+        Parameters
+        ----------
+        idx : int,
+            Band index. Indexation starts at 1
+
+        Returns
+        -------
+        min :
+            Band minimum
+        max :
+            Band maximum
+        mean :
+            Band mean
+        std :
+            Band deviation standard
+        """
+        ds = gdal.Open(self.filename, gdal.GA_ReadOnly)
+        band = ds.GetRasterBand(idx)
+
+        # False for approximation
+        return band.ComputeStatistics(False)
+
     def get_projection_from(self, raster):
         """TODO
         """
@@ -912,7 +971,76 @@ class Raster(Sized):
         else:
             return Raster(out_filename)
 
-    def rescale_bands(self, dstmin, dstmax, *idxs, **kw):
+    def rescale_bands(self, dstmin, dstmax, nodata_dst=None, out_filename=None,
+                      dtype=None,):
+        """
+        Rescales bands in the raster.
+
+        For each band, values are rescaled between given minimum and
+        maximum values.
+
+        Parameters
+        ----------
+        dstmin : float
+            Minimum value in the new scale.
+        dstmax : float
+            Maximum value in the new scale.
+        nodata_dst : float, (optional)
+            Nodata value to set in the output raster. The default nodata value
+            is used by default.
+        dtype: 'RasterDataType' (optionals),
+            Any data type supported by gdal, otb or numpy. By default,
+            'float32' is used.
+        out_filename : str
+            Path of the output file. If omitted, a default filename will be
+            chosen.
+
+        Returns
+        -------
+        `Raster` or None
+            Output raster or None if the raster is overwritten
+        """
+
+        out_filename = out_filename or add_suffix(self.filename,
+                                                  '_rescaled')
+        nodata_dst = nodata_dst or self.nodata_value
+        dtype = dtype or RasterDataType(lstr_dtype='float32')
+
+        # Apply rescaling on data different from nodata value
+        if self.nodata_value is not None:   # if there is nodata in raster
+            cond = '(im1b{i} == {nodata}) ? {nodata_dst} : ({formula})'
+        else:   # if there is no nodata value
+            cond = '{formula}'
+
+        expr = ''
+        for i in range(self.count):
+            band_min, band_max, _, _ = self.band_statistics(i + 1)
+            formula = ('{dstmin} ' +
+                       '+ (({dstmax} - {dstmin}) / ({srcmax} - {srcmin}))' +
+                       '* (im1b{i} - {srcmin})')
+            formula = formula.format(dstmin=dstmin, dstmax=dstmax,
+                                     srcmin=band_min, srcmax=band_max,
+                                     i=i + 1)
+            expr_b = cond.format(i=i+1, formula=formula,
+                                 nodata=self.nodata_value,
+                                 nodata_dst=nodata_dst)
+            if i + 1 < self.count:
+                expr_b += ';'
+            expr += expr_b
+
+        print 'BandMath expression is : ', expr
+
+        # Apply bandMathx application
+        BandMath = otb.Registry.CreateApplication("BandMathX")
+        BandMath.SetParameterStringList("il", [self.filename])
+        BandMath.SetParameterString("out", out_filename)
+        BandMath.SetParameterString("exp", expr)
+        BandMath.SetParameterOutputImagePixelType("out", dtype.otb_dtype)
+        BandMath.ExecuteAndWriteOutput()
+
+        return Raster(out_filename)
+
+    def rescale_bands_old(self, dstmin, dstmax, *idxs, **kw):
         """Rescales one or more bands in the raster.
 
         For each specified band, values are rescaled between given minimum and
@@ -968,6 +1096,69 @@ class Raster(Sized):
         else:
             return Raster(out_filename)
 
+    def normalize(self, nodata_dst=None, out_filename=None,
+                  dtype=None):
+        """
+        Normalize raster's band, i.e. remove the mean and divide by standard
+        deviation each band, without considering nodata values for statistics
+        computation.
+
+        Parameters
+        ----------
+        nodata_dst : float, (optional)
+            Nodata value to set in the output raster. The default nodata value
+            is used by default.
+        dtype: 'RasterDataType' (optionals),
+            Any data type supported by gdal, otb or numpy. By default,
+            'float32' is used.
+        out_filename : str
+            Path of the output file. If omitted, a default filename will be
+            chosen.
+
+        Returns
+        -------
+        `Raster` or None
+            Output raster
+        """
+
+        # Create an empty file with same size and dtype of float32
+        out_filename = out_filename or add_suffix(self.filename,
+                                                  '_normalized')
+        nodata_dst = nodata_dst or self.nodata_value
+        dtype = dtype or RasterDataType(lstr_dtype='float32')
+
+        # Apply rescaling on data different from nodata value if there is
+        # nodata
+        if self.nodata_value is not None:
+            cond = '(im1b{i} == {nodata}) ? {nodata_dst} : ({formula})'
+        else:
+            cond = '{formula}'
+
+        expr = ''
+        for i in range(self.count):
+            _, _, band_mean, band_std = self.band_statistics(i + 1)
+            formula = '(im1b{i} - {band_mean}) / {band_std}'
+            formula = formula.format(band_mean=band_mean, band_std=band_std,
+                                     i=i + 1)
+            expr_b = cond.format(i=i+1, formula=formula,
+                                 nodata=self.nodata_value,
+                                 nodata_dst=nodata_dst)
+            if i + 1 < self.count:
+                expr_b += ';'
+            expr += expr_b
+
+        print 'BandMath expression is : ', expr
+
+        # Apply bandMathx application
+        BandMath = otb.Registry.CreateApplication("BandMathX")
+        BandMath.SetParameterStringList("il", [self.filename])
+        BandMath.SetParameterString("out", out_filename)
+        BandMath.SetParameterString("exp", expr)
+        BandMath.SetParameterOutputImagePixelType("out", dtype.otb_dtype)
+        BandMath.ExecuteAndWriteOutput()
+
+        return Raster(out_filename)
+
     def fusion(self, pan, **kw):
         """
         Sharpen the raster with its corresponding panchromatic image.
@@ -1314,6 +1505,39 @@ class Raster(Sized):
 
         return out_raster
 
+    def row_col_to_lat_long(self, rows, cols):
+        """
+        Returns list x y coordinates from rows and cols position
+
+        Paremeters
+        ----------
+        rows : array like,
+            List of rows index
+        cols : array like,
+            List of columns index
+
+        Returns
+        -------
+        list_long : list of float
+            long coordinates
+        list_lat : list of float
+            lat coordinates
+        """
+
+        origin_x = self.transform[0]
+        origin_y = self.transform[3]
+        pixel_width = self.transform[1]
+        pixel_height = self.transform[5]
+
+        list_long = []
+        list_lat = []
+
+        for r, c in zip(rows, cols):
+            list_long.append(origin_x + pixel_width * c)
+            list_lat.append(origin_y + pixel_height * r)
+
+        return list_long, list_lat
+
     def _lsms_smoothing(self, spatialr, ranger, thres=0.1, rangeramp=0,
                         maxiter=10, **kw):
         """First step in a Large-Scale Mean-Shift (LSMS) segmentation: smooths