computeFeatures.py

import ConfigParser
import getopt
import glob
import json
import os
import platform
import subprocess
import sys
import warnings
import gdal
import xml.etree.ElementTree as ET
import otbApplication as otb

from mtdUtils import setNoDataValue


def readConfigFile(cfg_fn, section):
    config = ConfigParser.ConfigParser()
    config.read(cfg_fn)
    out = {}
    options = config.options(section)
    for option in options:
        out[option] = config.get(section, option)

    return out


def L8_readTile(fld):
    """
    mdf = glob.glob(fld + '/*_MTL*.txt')[0]
    loc_id = os.path.basename(os.path.split(mdf)[0])
    if 'MOSAIC' in loc_id:
        return 'MOSAIC'
    else:
        out = 'N.A.'
        with open(mdf) as mtl:
            for line in mtl:
                val = line.strip().split(' = ')
                if 'LANDSAT_SCENE_ID' in val[0]:
                    out = val[1][4:10]
        return out
    """
    # Rollback to folder name solution
    loc_id = os.path.basename(fld)
    out = loc_id.split('_')[-1]
    return out


def L8_readDate(fld):
    """
    mdf = glob.glob(fld + '/*_MTL*.txt')[0]
    loc_id = os.path.basename(os.path.split(mdf)[0])
    if 'MOSAIC' in loc_id:
        return loc_id[-15:-7]
    else:
        out = 'N.A.'
        with open(mdf) as mtl:
            for line in mtl:
                val = line.strip().split(' = ')
                if 'FILE_DATE' in val[0]:
                    out = val[1][0:10].replace('-', '')
        return out
    """
    # Rollback to folder name solution
    loc_id = os.path.basename(fld)
    out = loc_id.split('_')[-2]
    return out


def S2_readDate(fld):
    '''
    tli = glob.glob(fld + '/tileInfo.json')[0]
    loc_id = os.path.basename(os.path.split(tli)[0])
    if 'MOSAIC' in loc_id:
        return loc_id[-15:-7]
    else:
        with open(tli) as json_data:
            data = json.load(json_data)
        out = data['timestamp'][0:10].replace('-', '')
        return out
    '''
    # Rollback to folder name solution
    loc_id = os.path.basename(fld)
    out = loc_id.split('_')[-2]
    return out


def S2_readTile(fld):
    '''
    tli = glob.glob(fld + '/tileInfo.json')[0]
    loc_id = os.path.basename(os.path.split(tli)[0])
    if 'MOSAIC' in loc_id:
        out = 'MOSAIC'
    else:
        with open(tli) as json_data:
            data = json.load(json_data)
        out = str(data['utmZone']) + str(data['latitudeBand']) + str(data['gridSquare'])
    return out
    '''
    # Rollback to folder name solution
    loc_id = os.path.basename(fld)
    out = loc_id.split('_')[-1]
    return out


def S2THEIA_readDate(fld):
    '''
    mtd = glob.glob(fld + '/*_MTD_ALL.xml')[0]
    loc_id = os.path.basename(os.path.split(mtd)[0])
    if 'MOSAIC' in loc_id:
        return loc_id[-15:-7]
    else:
        dte = ET.parse(mtd).getroot().iter('ACQUISITION_DATE').next().text
        out = dte[0:4] + dte[5:7] + dte[8:10]
        return out
    '''
    # Rollback to folder name solution
    loc_id = os.path.basename(fld)
    out = loc_id.split('_')[-2]
    return out

def S2THEIA_readTile(fld):
    '''
    mtd = glob.glob(fld + '/*_MTD_ALL.xml')[0]
    loc_id = os.path.basename(os.path.split(mtd)[0])
    if 'MOSAIC' in loc_id:
        out = 'MOSAIC'
    else:
        out = ET.parse(mtd).getroot().iter('GEOGRAPHICAL_ZONE').next().text[1:]
    return out
    '''
    # Rollback to folder name solution
    loc_id = os.path.basename(fld)
    out = loc_id.split('_')[-1]
    return out


def VenusTHEIA_readDate(fld):
    '''
    mtd = glob.glob(fld + '/*_MTD_ALL.xml')[0]
    loc_id = os.path.basename(os.path.split(mtd)[0])
    if 'MOSAIC' in loc_id:
        return loc_id[-15:-7]
    else:
        dte = ET.parse(mtd).getroot().iter('ACQUISITION_DATE').next().text
        out = dte[0:4] + dte[5:7] + dte[8:10]
        return out
    '''
    # Rollback to folder name solution
    loc_id = os.path.basename(fld)
    out = loc_id.split('_')[-2]
    return out

def VenusTHEIA_readTile(fld):
    '''
    mtd = glob.glob(fld + '/*_MTD_ALL.xml')[0]
    loc_id = os.path.basename(os.path.split(mtd)[0])
    if 'MOSAIC' in loc_id:
        out = 'MOSAIC'
    else:
        out = ET.parse(mtd).getroot().iter('GEOGRAPHICAL_ZONE').next().text[1:]
    return out
    '''
    # Rollback to folder name solution
    loc_id = os.path.basename(fld)
    out = loc_id.split('_')[-1]
    return out


def checkCoregistered(f, cfg):
    coreg = False
    # Check if a VHR-coregistered version exists
    refl = None
    _refl = glob.glob(f + '/' + cfg['reflectance'])
    for x in _refl:
        if 'COREG' in x:
            refl = x
            coreg = True
    if refl is None:
        refl = _refl[0]

    gmask_name = None
    _gmask_name = glob.glob(f + '/' + cfg['gapmask'])
    if coreg:
        for x in _gmask_name:
            if 'COREG' in x:
                gmask_name = x
    else:
        gmask_name = _gmask_name[0]

    mask_name = None
    _mask_name = glob.glob(f + '/' + cfg['validitymask'])
    if coreg:
        for x in _mask_name:
            if 'COREG' in x:
                mask_name = x
    else:
        mask_name = _mask_name[0]

    #if refl is None or gmask_name is None or mask_name is None:
    if refl is None or mask_name is None:
        sys.exit("Problems with data or masks (check if fmask is performed).")

    return refl, gmask_name, mask_name


def L8_getExpressions(feat, cfg):
    expr = [None] * len(feat)
    
    bnds = ['B1','B2','B3','B4','B5','B6','B7','B9','B10','B11']
    for b in bnds:
        if b in feat:
            if cfg[b.lower()].isdigit():
                expr[feat.index(b)] = 'im1b' + cfg[b.lower()]
            else:
                warnings.warn("No " + b + "available on Landsat 8.")

    if 'NDVI' in feat:
        if cfg['b5'].isdigit() and cfg['b4'].isdigit():
            expr[feat.index('NDVI')] = '(im1b' + cfg['b5'] + ' - im1b' + cfg['b4'] + ') / (im1b' + cfg[
                'b5'] + ' + im1b' + cfg['b4'] + ')'
        else:
            warnings.warn("No NDVI produced on Landsat 8, invalid or non present B5 or B4.")

    if 'NDWI' in feat:
        if cfg['b3'].isdigit() and cfg['b5'].isdigit():
            expr[feat.index('NDWI')] = '(im1b' + cfg['b3'] + ' - im1b' + cfg['b5'] + ') / (im1b' + cfg[
                'b3'] + ' + im1b' + cfg['b5'] + ')'
        else:
            warnings.warn("No NDWI produced on Landsat 8, invalid or non present B3 or B5.")

    if 'MNDWI' in feat:
        if cfg['b3'].isdigit() and cfg['b6'].isdigit():
            expr[feat.index('MNDWI')] = '(im1b' + cfg['b3'] + ' - im1b' + cfg['b6'] + ') / (im1b' + cfg[
                'b3'] + ' + im1b' + cfg['b6'] + ')'
        else:
            warnings.warn("No MNDWI produced on Landsat 8, invalid or non present B3 or B6.")

    if 'SWNDVI' in feat:
        if cfg['b5'].isdigit() and cfg['b6'].isdigit():
            expr[feat.index('SWNDVI')] = '(im1b' + cfg['b5'] + ' - im1b' + cfg['b6'] + ') / (im1b' + cfg[
                'b5'] + ' + im1b' + cfg['b6'] + ')'
        else:
            warnings.warn("No SWNDVI produced on Landsat 8, invalid or non present B5 or B6.")

    if 'BRI' in feat:
        if cfg['b2'].isdigit() and cfg['b3'].isdigit() and cfg['b4'].isdigit() and cfg['b5'].isdigit() and cfg['b6'].isdigit() and cfg['b7'].isdigit():
            expr[feat.index('BRI')] = 'sqrt( im1b' + cfg['b2'] + ' * im1b' + cfg['b2'] + ' + im1b' + cfg[
                'b3'] + ' * im1b' + cfg['b3'] + ' + im1b' + cfg['b4'] + ' * im1b' + cfg['b4'] + ' + im1b' + cfg[
                                          'b5'] + ' * im1b' + cfg['b5'] + ' + im1b' + cfg['b6'] + ' * im1b' + cfg[
                                          'b6'] + ' + im1b' + cfg['b7'] + ' * im1b' + cfg['b7'] + ' )'
        else:
            warnings.warn("No BRI produced on Landsat 8, at least one invalid or non present band between B2-7.")

    return expr

def S2_getExpressions(feat, cfg):
    expr = [None] * len(feat)

    bnds = ['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B8A', 'B09', 'B10', 'B11', 'B12']
    for b in bnds:
        if b in feat:
            if cfg[b.lower()].isdigit():
                expr[feat.index(b)] = 'im1b' + cfg[b.lower()]
            else:
                warnings.warn("No " + b + "available on Sentinel-2.")

    if 'NDVI' in feat:
        if cfg['b08'].isdigit() and cfg['b04'].isdigit():
            expr[feat.index('NDVI')] = '(im1b' + cfg['b08'] + ' - im1b' + cfg['b04'] + ') / (im1b' + cfg[
                'b08'] + ' + im1b' + cfg['b04'] + ')'
        else:
            warnings.warn("No NDVI produced on Sentinel-2, invalid or non present B08 or B04.")

    if 'NDWI' in feat:
        if cfg['b03'].isdigit() and cfg['b08'].isdigit():
            expr[feat.index('NDWI')] = '(im1b' + cfg['b03'] + ' - im1b' + cfg['b08'] + ') / (im1b' + cfg[
                'b03'] + ' + im1b' + cfg['b08'] + ')'
        else:
            warnings.warn("No NDWI produced on Sentinel-2, invalid or non present B03 or B08.")

    if 'MNDWI' in feat:
        if cfg['b03'].isdigit() and cfg['b11'].isdigit():
            expr[feat.index('MNDWI')] = '(im1b' + cfg['b03'] + ' - im1b' + cfg['b11'] + ') / (im1b' + cfg[
                'b03'] + ' + im1b' + cfg['b11'] + ')'
        else:
            warnings.warn("No MNDWI produced on Sentinel-2, invalid or non present B03 or B11.")

    if 'SWNDVI' in feat:
        if cfg['b08'].isdigit() and cfg['b11'].isdigit():
            expr[feat.index('SWNDVI')] = '(im1b' + cfg['b08'] + ' - im1b' + cfg['b11'] + ') / (im1b' + cfg[
                'b08'] + ' + im1b' + cfg['b11'] + ')'
        else:
            warnings.warn("No SWNDVI produced on Sentinel-2, invalid or non present B08 or B11.")

    if 'BRI' in feat:
        if cfg['b02'].isdigit() and cfg['b03'].isdigit() and cfg['b04'].isdigit() and cfg['b05'].isdigit() and cfg['b06'].isdigit() and cfg['b07'].isdigit() and cfg['b08'].isdigit() and cfg['b8a'].isdigit() and cfg['b11'].isdigit() and cfg['b12'].isdigit():
            expr[feat.index('BRI')] = 'sqrt( im1b' + cfg['b02'] + ' * im1b' + cfg['b02'] + ' + im1b' + cfg[
                'b03'] + ' * im1b' + cfg['b03'] + ' + im1b' + cfg['b04'] + ' * im1b' + cfg['b04'] + ' + im1b' + cfg[
                                          'b05'] + ' * im1b' + cfg['b05'] + ' + im1b' + cfg['b06'] + ' * im1b' + cfg[
                                          'b06'] + ' + im1b' + cfg['b07'] + ' * im1b' + cfg['b07'] + ' + im1b' + cfg['b08'] + ' * im1b' + cfg['b08'] + ' + im1b' + cfg[
                                          'b8a'] + ' * im1b' + cfg['b8a'] + ' + im1b' + cfg['b11'] + ' * im1b' + cfg[
                                          'b11'] + ' + im1b' + cfg['b12'] + ' * im1b' + cfg['b12'] + ' )'
        else:
            warnings.warn("No BRI produced on Sentinel-2, at least one invalid or non present band among B02-B8A, B11 and B12.")

    # RED EDGE INDICES

    # Red-edge NDVI
    if 'NDRE' in feat:
        if cfg['b08'].isdigit() and cfg['b05'].isdigit():
            expr[feat.index('NDRE')] = '(im1b' + cfg['b08'] + ' - im1b' + cfg['b05'] + ') / (im1b' + cfg[
                'b08'] + ' + im1b' + cfg['b05'] + ')'
        else:
            warnings.warn("No NDRE produced on Sentinel-2, invalid or non present B08 or B05.")

    return expr

def S2THEIA_getExpressions(feat, cfg):
    expr = [None] * len(feat)

    bnds = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12']
    for b in bnds:
        if b in feat:
            if cfg[b.lower()].isdigit():
                expr[feat.index(b)] = 'im1b' + cfg[b.lower()]
            else:
                warnings.warn("No " + b + "available on Sentinel-2.")

    if 'NDVI' in feat:
        if cfg['b8'].isdigit() and cfg['b4'].isdigit():
            expr[feat.index('NDVI')] = '(im1b' + cfg['b8'] + ' - im1b' + cfg['b4'] + ') / (im1b' + cfg[
                'b8'] + ' + im1b' + cfg['b4'] + ')'
        else:
            warnings.warn("No NDVI produced on Sentinel-2, invalid or non present B8 or B4.")

    if 'NDWI' in feat:
        if cfg['b3'].isdigit() and cfg['b8'].isdigit():
            expr[feat.index('NDWI')] = '(im1b' + cfg['b3'] + ' - im1b' + cfg['b8'] + ') / (im1b' + cfg[
                'b3'] + ' + im1b' + cfg['b8'] + ')'
        else:
            warnings.warn("No NDWI produced on Sentinel-2, invalid or non present B3 or B8.")

    if 'MNDWI' in feat:
        if cfg['b3'].isdigit() and cfg['b11'].isdigit():
            expr[feat.index('MNDWI')] = '(im1b' + cfg['b3'] + ' - im1b' + cfg['b11'] + ') / (im1b' + cfg[
                'b3'] + ' + im1b' + cfg['b11'] + ')'
        else:
            warnings.warn("No MNDWI produced on Sentinel-2, invalid or non present B3 or B11.")

    if 'SWNDVI' in feat:
        if cfg['b8'].isdigit() and cfg['b11'].isdigit():
            expr[feat.index('SWNDVI')] = '(im1b' + cfg['b8'] + ' - im1b' + cfg['b11'] + ') / (im1b' + cfg[
                'b8'] + ' + im1b' + cfg['b11'] + ')'
        else:
            warnings.warn("No SWNDVI produced on Sentinel-2, invalid or non present B8 or B11.")

    if 'BRI' in feat:
        if cfg['b2'].isdigit() and cfg['b3'].isdigit() and cfg['b4'].isdigit() and cfg['b5'].isdigit() and cfg['b6'].isdigit() and cfg['b7'].isdigit() and cfg['b8'].isdigit() and cfg['b8a'].isdigit() and cfg['b11'].isdigit() and cfg['b12'].isdigit():
            expr[feat.index('BRI')] = 'sqrt( im1b' + cfg['b2'] + ' * im1b' + cfg['b2'] + ' + im1b' + cfg[
                'b3'] + ' * im1b' + cfg['b3'] + ' + im1b' + cfg['b4'] + ' * im1b' + cfg['b4'] + ' + im1b' + cfg[
                                          'b5'] + ' * im1b' + cfg['b5'] + ' + im1b' + cfg['b6'] + ' * im1b' + cfg[
                                          'b6'] + ' + im1b' + cfg['b7'] + ' * im1b' + cfg['b7'] + ' + im1b' + cfg['b8'] + ' * im1b' + cfg['b8'] + ' + im1b' + cfg[
                                          'b8a'] + ' * im1b' + cfg['b8a'] + ' + im1b' + cfg['b11'] + ' * im1b' + cfg[
                                          'b11'] + ' + im1b' + cfg['b12'] + ' * im1b' + cfg['b12'] + ' )'
        else:
            warnings.warn("No BRI produced on Sentinel-2, at least one invalid or non present band among B2-B8A, B11 and B12.")

    # RED EDGE INDICES

    # Red-edge NDVI
    if 'NDRE' in feat:
        if cfg['b8'].isdigit() and cfg['b5'].isdigit():
            expr[feat.index('NDRE')] = '(im1b' + cfg['b8'] + ' - im1b' + cfg['b5'] + ') / (im1b' + cfg[
                'b8'] + ' + im1b' + cfg['b5'] + ')'
        else:
            warnings.warn("No NDRE produced on Sentinel-2, invalid or non present B8 or B5.")

    return expr

def VenusTHEIA_getExpressions(feat, cfg):
    expr = [None] * len(feat)

    bnds = ['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B09', 'B10', 'B11', 'B12']
    for b in bnds:
        if b in feat:
            if cfg[b.lower()].isdigit():
                expr[feat.index(b)] = 'im1b' + cfg[b.lower()]
            else:
                warnings.warn("No " + b + "available on Venus.")

    if 'NDVI' in feat:
        if cfg['b11'].isdigit() and cfg['b07'].isdigit():
            expr[feat.index('NDVI')] = '(im1b' + cfg['b11'] + ' - im1b' + cfg['b07'] + ') / (im1b' + cfg[
                'b11'] + ' + im1b' + cfg['b07'] + ')'
        else:
            warnings.warn("No NDVI produced on Venus, invalid or non present B11 or B07.")

    if 'NDWI' in feat:
        if cfg['b04'].isdigit() and cfg['b11'].isdigit():
            expr[feat.index('NDWI')] = '(im1b' + cfg['b04'] + ' - im1b' + cfg['b11'] + ') / (im1b' + cfg[
                'b04'] + ' + im1b' + cfg['b11'] + ')'
        else:
            warnings.warn("No NDWI produced on Venus, invalid or non present B04 or B11.")

    if 'BRI' in feat:
        if cfg['b01'].isdigit() and cfg['b02'].isdigit() and cfg['b03'].isdigit() and cfg['b04'].isdigit() and cfg['b06'].isdigit() and cfg['b07'].isdigit() and cfg['b08'].isdigit() and cfg['b09'].isdigit() and cfg['b10'].isdigit() and cfg['b11'].isdigit() and cfg['b12'].isdigit():
            if 'B05' in feat:
                expr[feat.index('BRI')] = 'sqrt( im1b' + cfg['b01'] + ' * im1b' + cfg['b01'] + ' + im1b' + cfg['b02'] + ' * im1b' + cfg[
                                                             'b02'] + ' + im1b' + cfg['b03'] + ' * im1b' + cfg['b03'] + ' + im1b' + cfg[
                                                             'b04'] + ' * im1b' + cfg['b04'] + ' + im1b' + cfg['b05'] + ' * im1b' + cfg[
                                                             'b05'] + ' + im1b' + cfg['b06'] + ' * im1b' + cfg['b06'] + ' + im1b' + cfg[
                                                             'b07'] + ' * im1b' + cfg['b07'] + ' + im1b' + cfg['b08'] + ' * im1b' + cfg[
                                                             'b08'] + ' + im1b' + cfg['b09'] + ' * im1b' + cfg['b09'] + ' + im1b' + cfg[
                                                             'b10'] + ' * im1b' + cfg['b10'] + ' + im1b' + cfg['b11'] + ' * im1b' + cfg[
                                                             'b11'] + ' + im1b' + cfg['b12'] + ' * im1b' + cfg['b12'] + ' )'
            else :
                expr[feat.index('BRI')] = 'sqrt( im1b' + cfg['b01'] + ' * im1b' + cfg['b01'] + ' + im1b' + cfg['b02'] + ' * im1b' + cfg[
                                                             'b02'] + ' + im1b' + cfg['b03'] + ' * im1b' + cfg['b03'] + ' + im1b' + cfg[
                                                             'b04'] + ' * im1b' + cfg['b04'] + ' + im1b' + cfg['b06'] + ' * im1b' + cfg[
                                                             'b06'] + ' + im1b' + cfg['b07'] + ' * im1b' + cfg['b07'] + ' + im1b' + cfg[
                                                             'b08'] + ' * im1b' + cfg['b08'] + ' + im1b' + cfg['b09'] + ' * im1b' + cfg[
                                                             'b09'] + ' + im1b' + cfg['b10'] + ' * im1b' + cfg['b10'] + ' + im1b' + cfg[
                                                             'b11'] + ' * im1b' + cfg['b11'] + ' + im1b' + cfg['b12'] + ' * im1b' + cfg[
                                                             'b12'] + ' )'
        else:
            warnings.warn("No BRI produced on Venus, at least one invalid or non present band among B01-B12.")

    # RED EDGE INDICES

    # Red-edge NDVI
    if 'NDRE' in feat:
        if cfg['b11'].isdigit() and cfg['b08'].isdigit():
            expr[feat.index('NDRE')] = '(im1b' + cfg['b11'] + ' - im1b' + cfg['b08'] + ') / (im1b' + cfg[
                'b11'] + ' + im1b' + cfg['b08'] + ')'
        else:
            warnings.warn("No NDRE produced on Venus, invalid or non present B11 or B08.")

    return expr

def featureComputation(fld, cfg, sensor, out_fld = None, extent = None):
    if sensor not in ['S2', 'L8', 'S2_THEIA', 'Venus_THEIA']:
        sys.exit('Sensor not supported: ' + sensor)

    expr = []

    if cfg['feat'] == '':
        return
    else:
        feat = cfg['feat'].split(',')

    if cfg['gapfillingmode'] == 'BASE':
        gapf = 1 #Not implemented
    elif cfg['gapfillingmode'] == 'FEATURE':
        gapf = 2
    else:
        gapf = 0

    sh = False
    if platform.system() == 'Windows':
        sh = True

    supported = {}
    supported['L8'] = ['B1','B2','B3','B4','B5','B6','B7','B9','B10','B11','NDVI', 'NDWI', 'MNDWI', 'SWNDVI', 'BRI']
    supported['S2'] = ['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B8A', 'B09', 'B10', 'B11', 'B12', 'NDVI', 'NDWI', 'MNDWI', 'SWNDVI', 'BRI', 'NDRE']
    supported['S2_THEIA'] = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12', 'NDVI', 'NDWI', 'MNDWI', 'SWNDVI', 'BRI', 'NDRE']
    supported['Venus_THEIA'] = ['B01','B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B09', 'B10', 'B11', 'B12', 'NDVI', 'NDWI', 'BRI', 'NDRE']
    for ft in feat:
        if ft not in supported[sensor]:
            sys.exit("Unsupported feature for " + sensor + " selected: " + ft)

    if sensor == 'S2':
        tiles = [S2_readTile(x) for x in fld]
        expr = S2_getExpressions(feat, cfg)
    elif sensor == 'L8':
        tiles = [L8_readTile(x) for x in fld]
        expr = L8_getExpressions(feat, cfg)
    elif sensor == 'S2_THEIA':
        tiles = [S2THEIA_readTile(x) for x in fld]
        expr = S2THEIA_getExpressions(feat, cfg)
    elif sensor == 'Venus_THEIA':
        tiles = [VenusTHEIA_readTile(x) for x in fld]
        expr = VenusTHEIA_getExpressions(feat, cfg)

    minviews = cfg['minviews']

    u_tiles = list(set(tiles))
    src_tiles = dict((x, []) for x in u_tiles)
    feat_tiles = dict((x, []) for x in u_tiles)
    mask_tiles = dict((x, []) for x in u_tiles)
    vmask_tiles = dict((x, []) for x in u_tiles)
    out_sits = dict((x, []) for x in u_tiles)
    out_masks = dict((x, []) for x in u_tiles)
    out_nviews= dict((x, []) for x in u_tiles)
    out_sits_gapf = dict((x, []) for x in u_tiles)

    expr_tmp = ['(im2b1 * ' + x + ') + ((1-im2b1)*(-9999))' for x in expr]
    expr = '{' + ';'.join(expr_tmp) + '}'

    if out_fld is None:
        out_fld = fld[0]
        if len(fld) > 1:
            out_fld = fld[0] + '/../' + sensor + '_FEAT/'
            if not os.path.exists(out_fld):
                os.mkdir(out_fld)
    else:
        if not os.path.exists(out_fld):
            os.mkdir(out_fld)
        out_fld += '/' + sensor + '_FEAT/'
        if not os.path.exists(out_fld):
            os.mkdir(out_fld)

    envelope_fn = out_fld + '/' + sensor + '_envelope.tif'

    for f, t in zip(fld, tiles):
        refl, gmask_name, mask_name = checkCoregistered(f, cfg)

        # Build envelope for time series (if extent passed)
        if not os.path.exists(envelope_fn) and extent is not None and os.path.exists(extent):
            cmd = ['otbcli_Rasterization', '-in', extent, '-im', refl, '-out', envelope_fn, 'uint8', '-mode.binary.foreground', '1']
            subprocess.call(cmd, shell=sh)

        out_name = out_fld + os.path.basename(f) + '_FEAT.tif'
        #out_name_tmp = out_fld + os.path.basename(f) + '_TMP.tif'

        src_tiles[t].append(refl)
        feat_tiles[t].append(out_name)
        if os.path.exists(envelope_fn):
            gmask_name = out_fld + os.path.basename(f) + '_GAPS.tif'
            cmd = ['otbcli_BandMathX', '-il', envelope_fn, mask_name, '-exp', '{im1b1 != 0 && im2b1 == 0}', '-out', gmask_name, 'uint8']
            subprocess.call(cmd, shell=sh)
        mask_tiles[t].append(gmask_name)
        vmask_tiles[t].append(mask_name)
        #cmd = ['otbcli_BandMathX', '-il', refl, '-out', out_name_tmp, 'float', '-exp', expr]
        cmd = ['otbcli_BandMathX', '-il', refl, mask_name, '-out', out_name, 'float', '-exp', expr]
        subprocess.call(cmd, shell=sh)
        setNoDataValue(out_name,-9999)
        #cmd = ['otbcli_ManageNoData', '-in', out_name_tmp, '-out', out_name, 'float', '-mode', 'apply', '-mode.apply.mask', mask_name, '-mode.apply.ndval', '-9999']
        #subprocess.call(cmd, shell=sh)

        #os.remove(out_name_tmp)
        #if os.path.exists(out_name_tmp.replace('.tif', '.geom')):
        #    os.remove(out_name_tmp.replace('.tif', '.geom'))

    for l in feat_tiles:

        cmd = ['otbcli_ConcatenateImages', '-il'] + vmask_tiles[l] + ['-out',
                                                                     out_fld + sensor + '_' + l + '_VALMASK.tif',
                                                                     'uint8']
        subprocess.call(cmd, shell=sh)

        cmd = ['otbcli_ConcatenateImages', '-il'] + mask_tiles[l] + ['-out',
                                                                     out_fld + sensor + '_' + l + '_GAPMASK.tif',
                                                                     'uint8']
        subprocess.call(cmd, shell=sh)
        out_masks[l].append(out_fld + sensor + '_' + l + '_GAPMASK.tif')

        map_of_views = out_fld + sensor + '_' + l + '_NVIEWS.tif'
        ndates = len(mask_tiles[l])
        vexp = '(' + str(ndates) + ' - (' + ' + '.join(['im1b' + str(x) for x in range(1,ndates+1)]) + '))'
        cmd = ['otbcli_BandMathX', '-il', out_fld + sensor + '_' + l + '_GAPMASK.tif', '-exp', vexp, '-out', map_of_views, 'uint8']
        subprocess.call(cmd, shell=sh)
        out_nviews[l].append(map_of_views)

        for i in range(len(feat)):
            expr = '{' + ';'.join(['(-9999 * (1-im1b%d) + im%db%d * (im1b%d))' % (x + 1, x + 2, i + 1, x + 1) for x in range(len(feat_tiles[l]))]) + '}'

            # Shell version
            '''
            cmd = ['otbcli_BandMathX', '-il', out_fld + sensor + '_' + l + '_VALMASK.tif'] + feat_tiles[l] + ['-exp', expr, '-out',
                                                                 out_fld + sensor + '_' + l + '_' + feat[
                                                                     i] + '.tif', 'float']
            subprocess.call(cmd, shell=sh)
            '''
            # otpApplication version
            bmx = otb.Registry.CreateApplication('BandMathX')
            bmx.SetParameterStringList('il',[out_fld + sensor + '_' + l + '_VALMASK.tif'] + feat_tiles[l])
            bmx.SetParameterString('exp', expr)
            bmx.SetParameterString('out', out_fld + sensor + '_' + l + '_' + feat[i] + '.tif')
            bmx.ExecuteAndWriteOutput()

            setNoDataValue(out_fld + sensor + '_' + l + '_' + feat[i] + '.tif', -9999)
            out_sits[l].append(out_fld + sensor + '_' + l + '_' + feat[i] + '.tif')


            with open(out_fld + 'input_dates_' + l + '.txt', mode='w') as df:
                for f in src_tiles[l]:
                    ff = os.path.split(f)[0]
                    if sensor == 'L8':
                        df.write(L8_readDate(ff) + '\n')
                    elif sensor == 'S2':
                        df.write(S2_readDate(ff) + '\n')
                    elif sensor == 'S2_THEIA':
                        df.write(S2THEIA_readDate(ff) + '\n')
                    elif sensor == 'Venus_THEIA':
                        df.write(VenusTHEIA_readDate(ff) + '\n')
                        #for i in range(len(feat_tiles[l])):
                        #    os.remove(feat_tiles[l][i])
                        #    os.remove(feat_tiles[l][i].replace('.tif','.geom'))
    # Gapfilling
    if gapf == 2:
        for l in feat_tiles:
            ndates = len(src_tiles[l])
            for fn in out_sits[l]:
                ds = gdal.Open(fn)
                ncomp = ds.RasterCount / ndates
                ds = None
                cmd = ['otbcli_ImageTimeSeriesGapFilling', '-in', fn, '-mask', out_masks[l][0], '-comp', str(ncomp), '-it', 'linear', '-id', out_fld + 'input_dates_' + l + '.txt', '-out', fn.replace('.tif','_GAPF_tmp.tif'), 'float']
                subprocess.call(cmd,shell=sh)

                # Replace point with no views with no data
                cmd = ['otbcli_BandMath', '-il', out_nviews[l][0], '-exp', 'im1b1 >= ' + minviews, '-out', out_fld + '/tmpmsk.tif', 'uint8']
                subprocess.call(cmd, shell=sh)
                cmd = ['otbcli_ManageNoData', '-in', fn.replace('.tif','_GAPF_tmp.tif'), '-out', fn.replace('.tif','_GAPF.tif'), 'float', '-mode', 'apply', '-mode.apply.mask', out_fld + '/tmpmsk.tif', '-mode.apply.ndval', '-9999']
                subprocess.call(cmd, shell=sh)
                os.remove(out_fld + '/tmpmsk.tif')
                os.remove(fn.replace('.tif','_GAPF_tmp.tif'))
                if os.path.exists(fn.replace('.tif', '_GAPF_tmp.geom')):
                    os.remove(fn.replace('.tif', '_GAPF_tmp.geom'))

                out_sits_gapf[l].append(fn.replace('.tif','_GAPF.tif'))

    return out_sits, out_masks, out_sits_gapf, len(feat)


def main(argv):
    try:
        opts, args = getopt.getopt(argv, 'so:e:', ['series','output-folder=','extent='])
    except getopt.GetoptError as err:
        sys.exit(err)

    cfg_fn = args[0]
    fld = args[1]

    series = False
    out_fld = None
    extent = None

    for opt, val in opts:
        if opt in ['-s', '--series']:
            series = True
        elif opt in ['-o', '--output-folder']:
            out_fld = val
        elif opt in ['-e', '--extent']:
            extent = val

    L8_cfg = readConfigFile(cfg_fn, 'LANDSAT 8 CONFIGURATION')
    S2_cfg = readConfigFile(cfg_fn, 'SENTINEL-2 CONFIGURATION')
    S2THEIA_cfg = readConfigFile(cfg_fn, 'SENTINEL-2 THEIA CONFIGURATION')
    VenusTHEIA_cfg = readConfigFile(cfg_fn, 'VENUS THEIA CONFIGURATION')

    if series:
        flds = sorted([fld + '/' + f for f in next(os.walk(fld))[1]])
    else:
        flds = [fld]

    L8_flds = []
    S2_flds = []
    S2THEIA_flds = []
    VenusTHEIA_flds = []

    for f in flds:
        # check if S2 or L8 (search for metadata file)
        mdfL8 = glob.glob(f + '/*_MTL*.txt')
        mdfS2 = glob.glob(f + '/tileInfo.json')
        mdfS2THEIA = glob.glob(f + '/*_MTD_ALL.xml')
        mdfVenusTHEIA = glob.glob(f + '/*L2VALD.HDR')
        if len(mdfL8) > 0:
            L8_flds.append(f)
        elif len(mdfS2) > 0:
            S2_flds.append(f)
        elif len(mdfS2THEIA) > 0:
            S2THEIA_flds.append(f)
        elif len(mdfVenusTHEIA) > 0:
            VenusTHEIA_flds.append(f)


    if len(L8_flds) == 0 and len(S2_flds) == 0 and len(S2THEIA_flds) == 0 and len(VenusTHEIA_flds) == 0:
        sys.exit("No supported Landsat 8, Sentinel-2 or Venus data in " + fld)

    if len(L8_flds) > 0:
        featureComputation(L8_flds, L8_cfg, 'L8', out_fld=out_fld, extent=extent)
    if len(S2_flds) > 0:
        featureComputation(S2_flds, S2_cfg, 'S2', out_fld=out_fld, extent=extent)
    if len(S2THEIA_flds) > 0:
        featureComputation(S2THEIA_flds, S2THEIA_cfg, 'S2_THEIA', out_fld=out_fld, extent=extent)
    if len(VenusTHEIA_flds) > 0:
        featureComputation(VenusTHEIA_flds, VenusTHEIA_cfg, 'Venus_THEIA', out_fld=out_fld, extent=extent)

    return 0


if __name__ == '__main__':
    if len(sys.argv) < 2:
        sys.exit(
            'Usage: python computeFeatures.py [--series] [--output-folder <output folder>] [--extent <extent for full gap-filling>] <config-file> <image-dir>')
    main(sys.argv[1:])