genS2DownloadScript.py

import datetime
import getopt
import os
import platform
import stat
import subprocess
import sys
import glob

import ogr
import osr

from getSentinel2 import readS2TileInfo
from mtdUtils import setNoDataValue,randomword


def findtiles(shp):
    s2tg = ogr.Open('./aux_data/s2_tiling_grid.shp')
    clip = ogr.Open(shp)
    s2tg_ly = s2tg.GetLayer(0)
    clip_ly = clip.GetLayer(0)

    tiles = []

    s2tg_srs = s2tg_ly.GetSpatialRef()
    clip_srs = clip_ly.GetSpatialRef()
    tr = osr.CoordinateTransformation(clip_srs, s2tg_srs)

    cover = None
    covl = [0]

    for fr in s2tg_ly:
        clip_ly.ResetReading()
        gr = fr.GetGeometryRef()
        for f in clip_ly:
            g = f.GetGeometryRef()
            g.Transform(tr)
            if g.Intersects(gr):
                intsect = g.Intersection(gr)
                if cover is None:
                    cover = intsect
                else:
                    cover = cover.Union(intsect)
                covl.append(cover.GetArea() / g.GetArea())
                tl = fr.GetFieldAsString('Name')
                if tl not in tiles:
                    tiles.append(tl)

    dovl = [covl[i] - covl[i - 1] for i in range(1, len(covl))]

    idx = sorted(range(len(dovl)), key=lambda k: dovl[k], reverse=True)
    tot = 0.0
    i = 0
    while i < len(dovl) and tot < 1.0:
        tot += dovl[idx[i]]
        i += 1

    final_tiles = [tiles[k] for k in idx[0:i]]

    return final_tiles


def main(argv):
    try:
        opts, args = getopt.getopt(argv, 'o:crm:', ['min-coverage=', 'max-clouds=', 'scan-footprint'])
    except getopt.GetoptError as err:
        sys.exit(err)

    mindc = 0.0
    maxcc = 100.0
    clipping = False
    scan_foot = False

    # tile
    if os.path.exists(args[0]):
        tiles = findtiles(args[0])
    else:
        tiles = args[0].split(':')

    # Parse options
    opt_str = ''
    for opt, val in opts:
        if opt == '-o':
            opt_str = opt_str + ' -o \"' + val + '\"'
        elif opt == '-c' and os.path.exists(args[0]):
            clipping = True
            opt_str = opt_str + ' -c \"' + args[0] + '\"'
        elif opt == '-r':
            opt_str = opt_str + ' -r'
        elif opt == '-m':
            opt_str = opt_str + ' -m ' + val
        elif opt == '--min-coverage':
            mindc = float(val)
        elif opt == '--max-clouds':
            maxcc = float(val)
        elif opt == '--scan-footprint':
            scan_foot = True

    if platform.system() == 'Windows':
        f = open('S2DownloadScript.bat', 'w')
        f.write("@echo off\n")
        polycmd = 'gdal_polygonize '
    else:
        f = open('S2DownloadScript.sh', 'w')
        polycmd = 'gdal_polygonize.py '

    g = open('S2SearchLog.txt', 'w')

    DEVN = open(os.devnull, 'w')

    for tile in tiles:
        mrd = tile[0:2]
        tl1 = tile[2]
        tl2 = tile[3:5]
        # date and date-range
        sd = datetime.date(int(args[1][0:4]), int(args[1][4:6]), int(args[1][6:8]))
        ed = datetime.date(int(args[2][0:4]), int(args[2][4:6]), int(args[2][6:8])) + datetime.timedelta(1)
        cd = sd
        valid_dates = []
        print 'Scanning date interval for tile ' + tile + '...'
        while cd != ed:
            yy = str(cd.year)
            mm = str(cd.month)
            dd = str(cd.day)
            page = '/'.join([mrd, tl1, tl2, yy, mm, dd, '0'])
            url = 'http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/tiles/' + page + '/tileInfo.json'
            ret = subprocess.call('wget -q ' + url, shell=True)
            if ret == 0:
                # Patch: read over all possible "sequences" of a same date
                good_seq = 0
                good_info = readS2TileInfo('tileInfo.json')
                os.remove('tileInfo.json')
                curr_seq = 1
                while ret == 0:
                    page = page[:-1] + str(curr_seq)
                    url = 'http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/tiles/' + page + '/tileInfo.json'
                    ret = subprocess.call('wget -q ' + url, shell=True)
                    if ret == 0:
                        info = readS2TileInfo('tileInfo.json')
                        os.remove('tileInfo.json')
                        if info['dataCoveragePercentage'] > good_info['dataCoveragePercentage']:
                            good_info = info
                            good_seq = curr_seq
                    curr_seq += 1
                info = good_info

                datestr = ('%04d%02d%02d' % (cd.year, cd.month, cd.day), good_seq)
                if info['dataCoveragePercentage'] >= mindc and info['cloudyPixelPercentage'] <= maxcc:
                    # Check image footprint versus shape
                    clip_ok = True
                    if clipping and scan_foot:
                        fn = randomword(16)
                        page = page[:-1] + str(good_seq)
                        url = 'http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/tiles/' + page + '/B01.jp2'
                        subprocess.call('wget -q ' + url, shell=True)
                        subprocess.call('otbcli_BandMath -il B01.jp2 -exp \"im1b1>0\" -out ' + fn + '.tif uint8',
                                        shell=True, stdout=DEVN)
                        setNoDataValue(fn + '.tif')
                        subprocess.call(polycmd + fn + '.tif -f \"ESRI Shapefile" ' + fn + '.shp',
                                        shell=True, stdout=DEVN)
                        ds_query = ogr.Open(args[0])
                        ds_footp = ogr.Open(fn + '.shp')
                        ly_query = ds_query.GetLayer(0)
                        ly_footp = ds_footp.GetLayer(0)
                        ly_query_srs = ly_query.GetSpatialRef()
                        ly_footp_srs = ly_footp.GetSpatialRef()
                        ct = osr.CoordinateTransformation(ly_footp_srs, ly_query_srs)
                        intsec = False
                        for fq in ly_query:
                            for ff in ly_footp:
                                gf = ff.GetGeometryRef()
                                gf.Transform(ct)
                                if fq.GetGeometryRef().Intersects(gf):
                                    intsec = True
                        ds_query = None
                        ds_footp = None
                        for fl in glob.glob(fn + '.*'):
                            os.remove(fl)
                        os.remove('B01.jp2')
                        clip_ok = intsec

                    if clip_ok:
                        g.write('Adding    date ' + datestr[0] + ' for tile ' + tile + ' with ' + str(info['dataCoveragePercentage']) + '% data coverage and ' + str(info['cloudyPixelPercentage']) + '% cloudy pixels\n')
                        valid_dates.append(datestr)
                    else:
                        g.write('Rejecting date ' + datestr[0] + ' for tile ' + tile + ' with ' + str(
                            info['dataCoveragePercentage']) + '% data coverage and ' + str(
                            info['cloudyPixelPercentage']) + '% cloudy pixels (valid footprint does not intersect query)\n')
                else:
                    g.write('Rejecting date ' + datestr[0] + ' for tile ' + tile + ' with ' + str(
                        info['dataCoveragePercentage']) + '% data coverage and ' + str(
                        info['cloudyPixelPercentage']) + '% cloudy pixels\n')
            cd += datetime.timedelta(1)

        if len(valid_dates) > 0:
            print "Products have been found!"
            for dt in valid_dates:
                f.write('python getSentinel2.py ' + opt_str + ' ' + '--sequence ' + str(dt[1]) + ' ' + dt[0] + ' ' + tile + '\n')
    f.close()
    g.close()
    DEVN.close()

    if platform.system() == 'Linux':
        st = os.stat('S2DownloadScript.sh')
        os.chmod('S2DownloadScript.sh', st.st_mode | stat.S_IEXEC)


if __name__ == '__main__':
    if len(sys.argv) < 2:
        sys.exit(
            'Usage: python genS2DownloadScript.py [-o <output-dir>] [-c] [-r] [-m <bufsize>] [--min-coverage <perc>] [--max-clouds <perc>] [--scan-footprint] <shapefile OR colon separated tile numbers> <start date> <end date>')
    main(sys.argv[1:])