diff --git a/.gitignore b/.gitignore
index c6c5a741fc7899972605e0fff68b0ad42f8d5673..63d627dd6d736309587738385f7ca037d16ad0ea 100755
--- a/.gitignore
+++ b/.gitignore
@@ -60,3 +60,9 @@ data/l8_20130707\.tif\.aux\.xml
 data/l8_20130831\.tif\.aux\.xml
 
 qgis-auth\.db
+
+ymraster/\.spyderproject
+
+\.spyderworkspace
+
+\.spyderproject
diff --git a/ymraster/driver_ext.py b/ymraster/driver_ext.py
index 699c7ca955a50b92fbb9ec5414f6599b72d51384..60e5000f94df3e3cd1f1566803fb44a6b6d1873b 100755
--- a/ymraster/driver_ext.py
+++ b/ymraster/driver_ext.py
@@ -18,7 +18,8 @@ class DriverExt(GenericDriverExt):
 
     drivername_map = {'.tif': 'GTiff',
                       '.h5': 'HDF5',
-                      '.ers': 'ERS'}
+                      '.ers': 'ERS',
+                      '.vrt': 'VRT'}
 
     __slots__ = ()
 
diff --git a/ymraster/raster_dtype.py b/ymraster/raster_dtype.py
index 5caa6face7f6f02aa281e0a7a78c8f2109b945a1..a0c98d064dd2396dfc7fec1040bb1b5d3a737470 100755
--- a/ymraster/raster_dtype.py
+++ b/ymraster/raster_dtype.py
@@ -2,13 +2,22 @@
 
 try:
     import otbApplication as otb
+    OTB_IS_LOAD = True
 except ImportError as e:
-    raise ImportError(
-        str(e)
-        + "\n\nPlease install Orfeo Toolbox if it isn't installed yet.\n\n"
+    OTB_IS_LOAD = False
+    ERROR_OTB = ImportError(str(e)
+        + "\n\nThis functionnality is not available without Orfeo Toolbox. " +
+        "Please install Orfeo Toolbox if it isn't installed yet.\n\n"
         "Also, add the otbApplication module path "
         "(usually something like '/usr/lib/otb/python') "
-        "to the environment variable PYTHONPATH.\n")
+        "to the PYTHONPATH environment variable.")
+#    raise ImportError(
+#        str(e)
+#        + "\n\nPlease install Orfeo Toolbox if it isn't installed yet.\n\n"
+#        "Also, add the otbApplication module path "
+#        "(usually something like '/usr/lib/otb/python') "
+#        "to the environment variable PYTHONPATH.\n")
+
 try:
     from osgeo import gdal
 except ImportError as e:
@@ -32,81 +41,130 @@ class IdDefaultDict(dict):
         return self[key]
 
 
+#class DataType(object):
+#    """Abstract class for a data type (int16, int32, float32, etc.)"""
+#
+#    data_type_match = IdDefaultDict()
+#
+#    def __set__(self, instance, value):
+#        instance.otb_dtype = self.data_type_match[value]
+#
+#    def __get__(self, instance, owner):
+#        revert_match = {v: k for k, v in self.data_type_match.iteritems()}
+#        return revert_match[instance.otb_dtype]
+
 class DataType(object):
     """Abstract class for a data type (int16, int32, float32, etc.)"""
 
     data_type_match = IdDefaultDict()
 
     def __set__(self, instance, value):
-        instance.otb_dtype = self.data_type_match[value]
+        instance.gdal_dtype = self.data_type_match[value]
 
     def __get__(self, instance, owner):
         revert_match = {v: k for k, v in self.data_type_match.iteritems()}
-        return revert_match[instance.otb_dtype]
+        return revert_match[instance.gdal_dtype]
 
 
 class LStrDataType(DataType):
     """Represent a data type given in lower string format (eg. 'int16', 'int32',
     'float32', etc.)"""
 
-    data_type_match = {'uint8': otb.ImagePixelType_uint8,
-                       'uint16': otb.ImagePixelType_uint16,
-                       'uint32': otb.ImagePixelType_uint32,
-                       'int16': otb.ImagePixelType_int16,
-                       'int32': otb.ImagePixelType_int32,
-                       'float32': otb.ImagePixelType_float,
-                       'float64': otb.ImagePixelType_double}
+#    data_type_match = {'uint8': otb.ImagePixelType_uint8,
+#                       'uint16': otb.ImagePixelType_uint16,
+#                       'uint32': otb.ImagePixelType_uint32,
+#                       'int16': otb.ImagePixelType_int16,
+#                       'int32': otb.ImagePixelType_int32,
+#                       'float32': otb.ImagePixelType_float,
+#                       'float64': otb.ImagePixelType_double}
+
+    data_type_match = {'uint8': gdal.GDT_Byte,
+                       'uint16': gdal.GDT_UInt16,
+                       'uint32': gdal.GDT_UInt32,
+                       'int16': gdal.GDT_Int16,
+                       'int32': gdal.GDT_Int32,
+                       'float32': gdal.GDT_Float32,
+                       'float64': gdal.GDT_Float64}
 
 
 class UStrDataType(DataType):
     """Represent a data type given in upper string format (eg. 'Int16', 'Int32',
     'Float32', etc.)"""
 
-    data_type_match = {'UInt8': otb.ImagePixelType_uint8,
-                       'UInt16': otb.ImagePixelType_uint16,
-                       'UInt32': otb.ImagePixelType_uint32,
-                       'Int16': otb.ImagePixelType_int16,
-                       'Int32': otb.ImagePixelType_int32,
-                       'Float32': otb.ImagePixelType_float,
-                       'Float64': otb.ImagePixelType_double}
+#    data_type_match = {'UInt8': otb.ImagePixelType_uint8,
+#                       'UInt16': otb.ImagePixelType_uint16,
+#                       'UInt32': otb.ImagePixelType_uint32,
+#                       'Int16': otb.ImagePixelType_int16,
+#                       'Int32': otb.ImagePixelType_int32,
+#                       'Float32': otb.ImagePixelType_float,
+#                       'Float64': otb.ImagePixelType_double}
+    data_type_match = {'UInt8': gdal.GDT_Byte,
+                       'UInt16': gdal.GDT_UInt16,
+                       'UInt32': gdal.GDT_UInt32,
+                       'Int16': gdal.GDT_Int16,
+                       'Int32': gdal.GDT_Int32,
+                       'Float32': gdal.GDT_Float32,
+                       'Float64': gdal.GDT_Float64}
 
 
 class NumpyDataType(DataType):
     """Represent a data type for Numpy (eg. np.int16, np.int32, np.float32,
     etc.)"""
 
-    data_type_match = {np.uint8: otb.ImagePixelType_uint8,
-                       np.uint16: otb.ImagePixelType_uint16,
-                       np.uint32: otb.ImagePixelType_uint32,
-                       np.int16: otb.ImagePixelType_int16,
-                       np.int32: otb.ImagePixelType_int32,
-                       np.int64: otb.ImagePixelType_int32,
-                       np.float32: otb.ImagePixelType_float,
-                       np.float64: otb.ImagePixelType_double}
-
+#    data_type_match = {np.uint8: otb.ImagePixelType_uint8,
+#                       np.uint16: otb.ImagePixelType_uint16,
+#                       np.uint32: otb.ImagePixelType_uint32,
+#                       np.int16: otb.ImagePixelType_int16,
+#                       np.int32: otb.ImagePixelType_int32,
+#                       np.int64: otb.ImagePixelType_int32,
+#                       np.float32: otb.ImagePixelType_float,
+#                       np.float64: otb.ImagePixelType_double}
+    data_type_match = {np.uint8: gdal.GDT_Byte,
+                       np.uint16: gdal.GDT_UInt16,
+                       np.uint32: gdal.GDT_UInt32,
+                       np.int16: gdal.GDT_Int16,
+                       np.int32: gdal.GDT_Int32,
+                       np.float32: gdal.GDT_Float32,
+                       np.float64: gdal.GDT_Float64}
 
 class GdalDataType(DataType):
     """Represent a data type for gdal (eg. gdal.GDT_Int16, gdal.GDT_Iint32,
     gdal.GDT_Float32, etc.)"""
 
-    data_type_match = {gdal.GDT_Byte: otb.ImagePixelType_uint8,
-                       gdal.GDT_UInt16: otb.ImagePixelType_uint16,
-                       gdal.GDT_UInt32: otb.ImagePixelType_uint32,
-                       gdal.GDT_Int16: otb.ImagePixelType_int16,
-                       gdal.GDT_Int32: otb.ImagePixelType_int32,
-                       gdal.GDT_Float32: otb.ImagePixelType_float,
-                       gdal.GDT_Float64: otb.ImagePixelType_double}
+#    data_type_match = {gdal.GDT_Byte: otb.ImagePixelType_uint8,
+#                       gdal.GDT_UInt16: otb.ImagePixelType_uint16,
+#                       gdal.GDT_UInt32: otb.ImagePixelType_uint32,
+#                       gdal.GDT_Int16: otb.ImagePixelType_int16,
+#                       gdal.GDT_Int32: otb.ImagePixelType_int32,
+#                       gdal.GDT_Float32: otb.ImagePixelType_float,
+#                       gdal.GDT_Float64: otb.ImagePixelType_double}
+    def __set__(self, instance, value):
+        instance._gdal_type = value
+
+    def __get__(self, instance, owner):
+        return instance._gdal_type
 
 
 class OtbDataType(DataType):
     """Represent a data type for orfeo-toolbox
     (eg. otb.ImagePixelType_int16)"""
 
-    def __set__(self, instance, value):
-        instance._otb_type = value
-        
-    def __get__(self, instance, owner):
-        return instance._otb_type
+#    def __set__(self, instance, value):
+#        instance._otb_type = value
+#
+#    def __get__(self, instance, owner):
+#        return instance._otb_type
+    if OTB_IS_LOAD:
+        data_type_match = {otb.ImagePixelType_uint8: gdal.GDT_Byte,
+                           otb.ImagePixelType_uint16: gdal.GDT_UInt16,
+                           otb.ImagePixelType_uint32: gdal.GDT_UInt32,
+                           otb.ImagePixelType_int16: gdal.GDT_Int16,
+                           otb.ImagePixelType_int32: gdal.GDT_Int32,
+                           otb.ImagePixelType_float: gdal.GDT_Float32,
+                           otb.ImagePixelType_double: gdal.GDT_Float64}
+    else:
+        def __get__(self, instance, owner):
+            raise ERROR_OTB
 
 
 class RasterDataType(object):
@@ -124,7 +182,7 @@ class RasterDataType(object):
                  numpy_dtype=None,
                  otb_dtype=None,
                  gdal_dtype=None):
-                     
+
         if lstr_dtype:
             self.lstr_dtype = lstr_dtype
         elif ustr_dtype:
diff --git a/ymraster/ymraster.py b/ymraster/ymraster.py
index 7ca3744d67391dc249fa7505a3d6deceb15e4d1b..0114d7f72678771ea2644d91aceee35d3dcf897e 100755
--- a/ymraster/ymraster.py
+++ b/ymraster/ymraster.py
@@ -6,22 +6,29 @@ raster images.
 
 try:
     import otbApplication as otb
+    OTB_IS_LOAD = True
+    try:
+        app = otb.Registry.CreateApplication('Smoothing')
+        app.SetParameterString('out', 'foo.tif')
+    except AttributeError:
+        OTB_IS_LOAD = False
+        ERROR_OTB = AttributeError(
+            "This functionnality is not available without Orfeo Toolbox : "
+            "Unable to create otbApplication objects\n\n"
+            "Please set the ITK_AUTOLOAD_PATH environment variable "
+            "to the Orfeo Toolbox applications folder path "
+            "(usually something like '/usr/lib/otb/applications') ")
+        #raise ImportError(MSG_ERROR_OTB)
 except ImportError as e:
-    raise ImportError(
-        str(e)
-        + "\n\nPlease install Orfeo Toolbox if it isn't installed yet.\n\n"
+    OTB_IS_LOAD = False
+    ERROR_OTB = ImportError(str(e)
+        + "\n\nThis functionnality is not available without Orfeo Toolbox. " +
+        "Please install Orfeo Toolbox if it isn't installed yet.\n\n"
         "Also, add the otbApplication module path "
         "(usually something like '/usr/lib/otb/python') "
         "to the PYTHONPATH environment variable.")
-try:
-    app = otb.Registry.CreateApplication('Smoothing')
-    app.SetParameterString('out', 'foo.tif')
-except AttributeError:
-    raise ImportError(
-        "Unable to create otbApplication objects\n\n"
-        "Please set the ITK_AUTOLOAD_PATH environment variable "
-        "to the Orfeo Toolbox applications folder path "
-        "(usually something like '/usr/lib/otb/applications') ")
+    #raise ImportError(MSG_ERROR)
+
 
 try:
     from osgeo import osr, gdal
@@ -49,6 +56,16 @@ from time import mktime
 import os
 import shutil
 from tempfile import gettempdir
+from functools import wraps
+
+def otb_function(f):
+    @wraps(f)
+    def wrapper(*args, **kwds):
+        if not OTB_IS_LOAD:
+            raise ERROR_OTB
+        return f(*args, **kwds)
+    return wrapper
+
 
 def add_suffix(path_file, suffix):
     """
@@ -95,53 +112,59 @@ def write_file(out_filename, array=None, overwrite=True,
     """Writes a NumPy array to an image file.
 
     If there is no array (array is None), the function simply create an empty
-    image file of given size (width, height, depth). In other words, if array is
-    None, width, height and depth must be specified.
+    image file of given size (width, height, depth). In other words, if array
+    is None, width, height and depth must be specified.
 
     If array is not None, then all other parameters are optional: if the file
-    exists, array will be written into the file. If the file does not exists, a
-    new file will be created with same size and type than the array.
+    exists, array will be written into the file. If the file does not exists,
+    a new file will be created with same size and type than the array.
 
     This allows to write a file block by block. First create an empty file of
-    correct size and type. Then write each block into the file using the xoffset
-    and yoffset parameters.
+    correct size and type. Then write each block into the file using the
+    xoffset and yoffset parameters.
 
-    :param out_filename: path to the output file
-    :type out_filename: str
-    :param array: the NumPy array to save. None means that an empty file will be
-                  created or, if file exists, that no data will be written
-                  (except metadata if specified)
-    :type array: np.ndarray
-    :param overwrite: if True, overwrite file if exists. If False, open the
-                    file, then update it or create it if it does not exist.
-                    True by default.
-    :type overwrite: bool
-    :param xoffset: horizontal offset. First index from which to write the array
-                    in the output file if the array is smaller (default: 0)
-    :type xoffset: float
-    :param yoffset: Vertical offset. First index from which to write the array
-                    in the output file if the array is smaller (default: 0)
-    :type yoffset: float
-    :param band_idx: depth offset. First band from which to write the array
-                    in the output file if the array is smaller (default: 1)
-    :type band_idx: int
-    :param width: horizontal size of the image to be created
-    :type width: int
-    :param height: vertical size of the image to be created
-    :type width: int
-    :param depth: number of bands of the image to be created
-    :type depth: int
-    :param dtype: data type to use for the output file. None means that the file
-                  already exists
-    :type dtype: RasterDataType
-    :param date_time: date/time to write in the output file metadata
-    :type date_time: datetime.datetime
-    :param nodata_value: value to set as NODATA in the output file metadata
-    :type nodata_value: dtype
-    :param srs: projection to write in the output file metadata
-    :type srs: osr.SpatialReference
-    :param transform: geo-transformation to use for the output file
-    :type transform: 6-tuple of floats
+    Parameters
+    ----------
+    out_filename : str,
+        Path to the output file.
+    array : np.ndarray,
+        The NumPy array to save. None means that an empty file will be
+        created or, if file exists, that no data will be written (except
+        metadata if specified)
+    overwrite : bool,
+        if True, overwrite file if exists. If False, open the
+        file, then update it or create it if it does not exist. True by
+        default.
+    xoffset : flaot,
+        Horizontal offset. First index from which to write the array
+        in the output file if the array is smaller (default: 0).
+    yoffset : float,
+        Vertical offset. First index from which to write the array
+        in the output file if the array is smaller (default: 0).
+    band_idx : int
+        Depth offset. First band from which to write the array
+        in the output file if the array is smaller (default: 1).
+    kw
+    --
+    width : int,
+        Horizontal size of the image to be created (in pixel)
+    height : int,
+        Vertical size of the image to be created (in pixel)
+    depth : int,
+        Number of bands of the image to be created.
+    dtype : RasterDataType
+        Data type to use for the output file. None means that the file
+        already exists.
+    date_time : datetime.datetime
+        Date/time to write in the output file metadata.
+    nodata_value : dtype
+        Value to set as NODATA in the output file metadata.
+    srs : osr.SpatialReference
+        projection to write in the output file metadata.
+    transform : 6-tuple of floats
+        Geo-transformation to use for the output file.
+    band_info : list of string:
+        Info of each band to write in the output file metadata
     """
     # Size & data type of output image
     xsize, ysize = (kw['width'], kw['height']) \
@@ -173,7 +196,8 @@ def write_file(out_filename, array=None, overwrite=True,
     # Set metadata
     if kw.get('date_time'):
         out_ds.SetMetadata(
-            {'TIFFTAG_DATETIME': kw['date_time'].strftime('%Y:%m:%d %H:%M:%S')})
+            {'TIFFTAG_DATETIME':
+             kw['date_time'].strftime('%Y:%m:%d %H:%M:%S')})
     if kw.get('srs'):
         out_ds.SetProjection(kw['srs'].ExportToWkt())
     if kw.get('transform'):
@@ -188,6 +212,10 @@ def write_file(out_filename, array=None, overwrite=True,
         band.WriteArray(array, xoff=xoffset, yoff=yoffset)
         if kw.get('nodata_value'):
             band.SetNoDataValue(kw.get('nodata_value'))
+        if kw.get('band_info'):
+            metadata = band.GetMetadata()
+            metadata['band_info'] = kw.get('band_info')[0]
+            band.SetMetadata(metadata)
         band.FlushCache()
     else:
         for i in range(number_bands):
@@ -195,11 +223,15 @@ def write_file(out_filename, array=None, overwrite=True,
             band.WriteArray(array[:, :, i], xoff=xoffset, yoff=yoffset)
             if kw.get('nodata_value'):
                 band.SetNoDataValue(kw.get('nodata_value'))
+            if kw.get('band_info'):
+                metadata = band.GetMetadata()
+                metadata['band_info'] = kw.get('band_info')[i]
+                band.SetMetadata(metadata)
             band.FlushCache()
     band = None
     out_ds = None
 
-
+@otb_function
 def concatenate_rasters_old(*rasters, **kw):
     """Write a raster which is the concatenation of the given rasters, in order.
 
@@ -256,6 +288,8 @@ def concatenate_rasters_old(*rasters, **kw):
         shutil.copy(out_filename, raster0.filename)
         os.remove(out_filename)
 
+
+@otb_function
 def concatenate_rasters(rasters, **kw):
     """Write a raster which is the concatenation of the given rasters, in order.
 
@@ -637,6 +671,125 @@ class Raster(Sized):
         # Close file
         ds = None
 
+    def write_metadata(self, metadata):
+        """
+        Write metada in the file
+
+        Parameters
+        ----------
+        metadata : dictionnary
+            Key are Tag of the metadata, values is the corresponding
+            description
+
+        Limitations
+        -----------
+        Works for Tif files only.
+
+        """
+        ds = gdal.Open(self._filename, gdal.GA_Update)
+
+        metadata_orign = ds.GetMetadata()
+        metadata_orign.update(metadata)
+        ds.SetMetadata(metadata_orign)
+        ds = None
+        self.refresh()
+
+
+    def write_band_metadata(self, idx, metadata):
+        """
+        Write band metada in the file
+
+        Parameters
+        ----------
+        idx : int
+            Band number. Indexation starts at 1.
+        metadata : dictionnary
+            Key are Tag of the metadata, values is the corresponding
+            description
+
+        Limitations
+        -----------
+        Works for Tif files only.
+
+        """
+        ds = gdal.Open(self._filename, gdal.GA_Update)
+        band = ds.GetRasterBand(idx)
+        metadata_orign = band.GetMetadata()
+        metadata_orign.update(metadata)
+        band.SetMetadata(metadata)
+        ds = None
+        self.refresh()
+
+    def get_band_metadata(self, idx):
+        """
+        Get band metada in the file
+
+        Parameters
+        ----------
+        idx : int
+            Band number. Indexation starts at 1.
+
+        Returns
+        -------
+        metadata : dictionnary
+            Key are Tag of the metadata, values is the corresponding
+            description
+        """
+        ds = gdal.Open(self._filename, gdal.GA_Update)
+        band = ds.GetRasterBand(idx)
+        metadata = band.GetMetadata()
+        ds = None
+        return metadata
+
+    @otb_function
+    def change_dtype(self, out_dtype, out_filename):
+        """
+        Convert Raster datatype to an other using Otb BandMath application
+
+        Parameters
+        ----------
+        dytpe : str
+            Can be one of this type :
+                - 'uint8'
+                - 'uint16'
+                - 'uint32'
+                - 'int16'
+                - 'int32'
+                - 'float32'
+                - 'float64'
+        Returns
+        -------
+        out_rst : `Raster`
+            Output raster
+
+        Limitation
+        ----------
+        Have not been tested in case rescaling and.or resampling is needed
+        """
+
+        dtype = RasterDataType(lstr_dtype=out_dtype).otb_dtype
+        expr = ''
+        for i in range(self.count):
+            expr_b = 'im1b{}'.format(i + 1)
+            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)
+        BandMath.ExecuteAndWriteOutput()
+
+        # write new nodata value in metadata
+        out_rst = Raster(out_filename)
+
+        return out_rst
+
     def band_statistics(self, idx):
         """
         Compute band statistic, without taking into account nodata values
@@ -663,6 +816,72 @@ class Raster(Sized):
         # False for approximation
         return band.ComputeStatistics(False)
 
+    def compute_percentile(self, list_per, list_idx=None,):
+        """
+        Compute and return percentile(s) of Raster's bands and write it in
+        bands' metadata.
+
+        Parameters
+        ---------
+        list_per : list of int.
+            Percentiles must be comprised between 0 and 100
+        list_idx: list of int (optional)
+            List of band index (idexation starts at 1) on which compute
+            percentile(s). All band are processed by default.
+
+        Returns
+        -------
+        metadata_all_bands : dict of dict:
+            Dictionnary of metadatas for each band.
+        """
+        metadata_all_bands = {}
+        ds = gdal.Open(self.filename, gdal.GA_ReadOnly)
+
+        if list_idx is None:
+            list_idx = range(1, self.count + 1)
+
+        for idx in list_idx:
+            tab = self.array_from_bands(idx, mask_nodata=True)
+            tab = tab.compressed()
+            band = ds.GetRasterBand(idx)
+            metadata = band.GetMetadata()
+            for per in list_per:
+                per_label = 'STATISTICS_PER_{}'.format(per)
+                per_value = np.nanpercentile(tab, per)
+                metadata[per_label] = str(per_value)
+            metadata_all_bands[idx] = metadata
+            band.SetMetadata(metadata)
+
+        return metadata_all_bands
+
+    def get_percentile(self, per, idx):
+        """
+        Return percentile of a Raster's band according to its metadata.
+
+        Parameters
+        ---------
+        per : int.
+            Percentile, must be comprised between 0 and 100.
+        idx: int
+            Band index (idexation starts at 1) on which compute
+            percentile `per`.
+
+        Returns
+        -------
+        percentile_value : float
+            Percentile value
+
+        See also
+        --------
+        Raster.compute_percentile()
+        """
+
+        ds = gdal.Open(self.filename, gdal.GA_ReadOnly)
+        band = ds.GetRasterBand(idx)
+        per_label = 'STATISTICS_PER_{}'.format(per)
+
+        return float(band.GetMetadata()[per_label])
+
     def get_projection_from(self, raster):
         """TODO
         """
@@ -850,6 +1069,7 @@ class Raster(Sized):
                    block_win[0],
                    block_win[1])
 
+    @otb_function
     def remove_bands(self, *idxs, **kw):
         """
         Saves a new raster with the specified bands removed.
@@ -906,6 +1126,7 @@ class Raster(Sized):
         else:
             return Raster(out_filename)
 
+    @otb_function
     def extract_band(self,*idxs,**kw):
         """
         Saves a new raster with the specified bands extracted and the others
@@ -938,18 +1159,19 @@ class Raster(Sized):
         SplitImage.ExecuteAndWriteOutput()
 
         # Out file
-        out_filename =  kw.get('out_filename',os.path.join(gettempdir(), 'bands_removed.tif'))
+        out_filename = kw.get('out_filename',
+                              os.path.join(gettempdir(), 'bands_removed.tif'))
 
         # Concatenate the mono-band images with only the wanted band(s)
         #------------------------------------------------------------------
-        list_path = [os.path.join(gettempdir(), 'splitted_{}.tif'.format(i))
-                     for i in range(self._count)
-                     if i + 1 in indices]
+        list_path = [os.path.join(gettempdir(),
+                                  'splitted_{}.tif'.format(i - 1))
+                     for i in indices]
 
         #Re order if indicated
-        temp_list = [(idx, path) for idx, path in zip(indices,list_path)]
-        temp_list.sort()
-        list_path = [path for (idx, path) in temp_list]
+#        temp_list = [(idx, path) for idx, path in zip(indices,list_path)]
+#        temp_list.sort()
+#        list_path = [path for (idx, path) in temp_list]
 
         #concatenate here
         ConcatenateImages = otb.Registry.CreateApplication("ConcatenateImages")
@@ -971,6 +1193,7 @@ class Raster(Sized):
         else:
             return Raster(out_filename)
 
+    @otb_function
     def rescale_bands(self, dstmin, dstmax, nodata_dst=None, out_filename=None,
                       dtype=None,):
         """
@@ -1096,6 +1319,7 @@ class Raster(Sized):
         else:
             return Raster(out_filename)
 
+    @otb_function
     def normalize(self, nodata_dst=None, out_filename=None,
                   dtype=None):
         """
@@ -1157,8 +1381,14 @@ class Raster(Sized):
         BandMath.SetParameterOutputImagePixelType("out", dtype.otb_dtype)
         BandMath.ExecuteAndWriteOutput()
 
-        return Raster(out_filename)
+        # write new nodata value in metadata
+        out_rst = Raster(out_filename)
+        if nodata_dst:
+            out_rst.nodata_value = nodata_dst
+
+        return out_rst
 
+    @otb_function
     def fusion(self, pan, **kw):
         """
         Sharpen the raster with its corresponding panchromatic image.
@@ -1208,6 +1438,7 @@ class Raster(Sized):
         else:
             return Raster(out_filename)
 
+    @otb_function
     @fix_missing_proj
     def radiometric_indices(self, *indices, **kw):
         """Saves a raster of radiometric indices about the raster.
@@ -1281,6 +1512,7 @@ class Raster(Sized):
             out_raster.date_time = self._date_time
         return out_raster
 
+    @otb_function
     def ndvi(self, red_idx, nir_idx, **kw):
         """Saves the Normalized Difference Vegetation Index (NDVI) of the
         raster.
@@ -1308,6 +1540,7 @@ class Raster(Sized):
                                         nir_idx=nir_idx,
                                         out_filename=out_filename)
 
+    @otb_function
     def ndwi(self, nir_idx, mir_idx, **kw):
         """Saves the Normalized Difference Water Index (NDWI) of the raster.
 
@@ -1334,6 +1567,7 @@ class Raster(Sized):
                                         mir_idx=mir_idx,
                                         out_filename=out_filename)
 
+    @otb_function
     def mndwi(self, green_idx, mir_idx, **kw):
         """Saves the Modified Normalized Difference Water Index (MNDWI) of the
         raster.
@@ -1363,6 +1597,7 @@ class Raster(Sized):
 
     ndsi = mndwi
 
+    @otb_function
     def append(self, *rasters):
         """Append the given rasters into the current one.
 
@@ -1374,6 +1609,7 @@ class Raster(Sized):
         raster_list = [self] + list(rasters)
         concatenate_rasters(raster_list)
 
+    @otb_function
     def apply_mask(self, mask_raster, mask_value=1, set_value=None, **kw):
         """Apply a mask to the raster: set all masked pixels to a given value.
 
@@ -1442,7 +1678,7 @@ class Raster(Sized):
         masked_bands = [Raster(os.path.join(gettempdir(),
                                             'mono_masked_{}.tif'.format(i)))
                         for i in range(self._count)]
-        concatenate_rasters(*masked_bands, out_filename=out_filename)
+        concatenate_rasters(masked_bands, out_filename=out_filename)
         out_raster = Raster(out_filename)
         out_raster.nodata_value = set_value
         print set_value
@@ -1453,6 +1689,7 @@ class Raster(Sized):
 
         return out_raster
 
+    @otb_function
     def change_value(self, dict_change, idx_band, **kw):
         """
         Change values by others in raster band.
@@ -1538,6 +1775,7 @@ class Raster(Sized):
 
         return list_long, list_lat
 
+    @otb_function
     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
@@ -1599,6 +1837,7 @@ class Raster(Sized):
 
         return Raster(out_filename), Raster(out_spatial_filename)
 
+    @otb_function
     def _lsms_segmentation(self, spatialr, ranger, spatial_raster,
                            block_size=None, **kw):
         """Second step in a LSMS segmentation: performs the actual object
@@ -1669,6 +1908,7 @@ class Raster(Sized):
 
         return Raster(out_filename)
 
+    @otb_function
     @fix_missing_proj
     def _lsms_merging(self, object_minsize, smoothed_raster, block_size=None,
                       **kw):
@@ -1735,6 +1975,7 @@ class Raster(Sized):
         else:
             return Raster(out_filename)
 
+    @otb_function
     def _lsms_vectorization(self, orig_raster, block_size=None, **kw):
         """Fourth and Last (optional) step in a LSMS segmentation: vectorize a
         labeled segmented image, turn each object into a polygon. Each polygon
@@ -1785,6 +2026,7 @@ class Raster(Sized):
         LSMSVectorization.SetParameterInt("tilesizey", tilesizey)
         LSMSVectorization.ExecuteAndWriteOutput()
 
+    @otb_function
     def lsms_segmentation(self,
                           spatialr,
                           ranger,