computeFeatures.py

import ConfigParser
import getopt
import glob
import json
import os
import platform
import subprocess
import sys
import warnings
import gdal

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]

    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


def L8_readDate(fld):
    mdf = glob.glob(fld + '/*_MTL*.txt')[0]

    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


def S2_readDate(fld):
    tli = glob.glob(fld + '/tileInfo.json')[0]
    with open(tli) as json_data:
        data = json.load(json_data)

    out = data['timestamp'][0:10].replace('-', '')

    return out


def S2_readTile(fld):
    tli = glob.glob(fld + '/tileInfo.json')[0]
    with open(tli) as json_data:
        data = json.load(json_data)

    out = str(data['utmZone']) + str(data['latitudeBand']) + str(data['gridSquare'])

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

    if 'MNDWI' in feat:
        if cfg['b5'].isdigit() and cfg['b6'].isdigit():
            expr[feat.index('MNDWI')] = '(im1b' + cfg['b5'] + ' - im1b' + cfg['b6'] + ') / (im1b' + cfg[
                'b5'] + ' + im1b' + cfg['b6'] + ')'
        else:
            warnings.warn("No MNDWI 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 'MNDVI' in feat:
        if cfg['b03'].isdigit() and cfg['b11'].isdigit():
            expr[feat.index('MNDVI')] = '(im1b' + cfg['b03'] + ' - im1b' + cfg['b11'] + ') / (im1b' + cfg[
                'b03'] + ' + im1b' + cfg['b11'] + ')'
        else:
            warnings.warn("No MNDVI produced on Sentinel-2, invalid or non present B03 or B11.")

    if 'MNDWI' in feat:
        if cfg['b08'].isdigit() and cfg['b11'].isdigit():
            expr[feat.index('MNDWI')] = '(im1b' + cfg['b08'] + ' - im1b' + cfg['b11'] + ') / (im1b' + cfg[
                'b08'] + ' + im1b' + cfg['b11'] + ')'
        else:
            warnings.warn("No MNDWI 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 'RNDVI' in feat:
        if cfg['b08'].isdigit() and cfg['b06'].isdigit():
            expr[feat.index('RNDVI')] = '(im1b' + cfg['b08'] + ' - im1b' + cfg['b06'] + ') / (im1b' + cfg[
                'b08'] + ' + im1b' + cfg['b06'] + ')'
        else:
            warnings.warn("No RNDVI produced on Sentinel-2, invalid or non present B08 or B06.")

    return expr


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

    expr = []

    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', 'MNDVI', 'MNDWI', 'BRI']
    supported['S2'] = ['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B8A', 'B09', 'B10', 'B11', 'B12', 'NDVI', 'NDWI', 'MNDVI', 'MNDWI', 'BRI', 'RNDVI']
    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)

    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 = '{' + ';'.join(expr) + '}'

    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]
        #subprocess.call(cmd, shell=sh)
        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]))]) + '}'
            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)
            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')
                    #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'))

                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, 's', ['series'])
    except getopt.GetoptError as err:
        sys.exit(err)

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

    series = False

    for opt, val in opts:
        if opt in ['-s', '--series']:
            series = True

    L8_cfg = readConfigFile(cfg_fn, 'LANDSAT 8 CONFIGURATION')
    S2_cfg = readConfigFile(cfg_fn, 'SENTINEL-2 CONFIGURATION')

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

    L8_flds = []
    S2_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')
        if len(mdfL8) > 0:
            L8_flds.append(f)
        elif len(mdfS2) > 0:
            S2_flds.append(f)

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

    if len(L8_flds) > 0:
        featureComputation(L8_flds, L8_cfg, 'L8')
    if len(S2_flds) > 0:
        featureComputation(S2_flds, S2_cfg, 'S2')

    return 0


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