awdc.py

#!/g/data/v10/public/modules/agdc-py3-env/20171214/envs/agdc/bin/python

import pyproj
import argparse
import glob
import os
import sys
import shutil
import uuid
import subprocess
from datetime import datetime

def get_agdc_files(lon, lat):
    gda94aa = pyproj.Proj(init='epsg:3577')
    wgs84 = pyproj.Proj(init='epsg:4326')
    x, y = pyproj.transform(wgs84, gda94aa, lon, lat)

    print(x, y)
    print("/g/data/rs0/datacube/002/LS8_OLI_NBAR/{}_{}/*".format(int(x/100000), int(y/100000)))
    return sorted(glob.glob("/g/data/rs0/datacube/002/LS8_OLI_NBAR/{}_{}/*".format(int(x/100000), int(y/100000))))



if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="""Modis Vegetation Analysis argument parser""")
    parser.add_argument('-p', '--point', type=str, required=True, help="Specify the requested extent as [lon,lat].")
    parser.add_argument('-o', '--output', required=True, type=str, help="Full path to destination.")
    parser.add_argument('-tmp', '--tmp', required=True, type=str, help="Full path to destination.")
    args = parser.parse_args()

    print(args.point)
    lon, lat = [float(v) for v in args.point.split(' ')]

    band = "{}/band.nc".format(args.tmp)
    band_time = "{}/band_time.nc".format(args.tmp)
    merged = "{}/merged.nc".format(args.tmp)
    files = get_agdc_files(lon, lat)

    for n, f in enumerate(files):
        if n == 10:
            sys.exit(0)

        gdal_h = 'NETCDF:"{}":red'.format(f)
        tmp_file = os.path.join(uuid.uuid4().hex)
        os.system("gdalwarp -of netCDF -r cubic -t_srs '+proj=longlat +datum=WGS84 +no_defs' -te 149.0 -35.0 149.05 -34.95 -tr 0.01 -0.01 {} {}".format(gdal_h, tmp_file))

        result = subprocess.run(["cdo", "showdate", "{}".format(f)], stdout=subprocess.PIPE)
        dates = result.stdout.rstrip().decode("utf-8").split()
        for b, d in enumerate(dates):
            print(d, b)
            os.system("cdo -selvar,Band{} {} {}".format(b+1, tmp_file, band))
            os.system("cdo -setdate,{} {} {}".format(d, band, band_time))
            os.remove(band)
            os.system("cdo showdate {}".format(band_time))
            if not os.path.isfile(args.output):
                shutil.move(band_time, args.output)
                continue

            print("MERGE")
            os.system("cdo mergetime {} {} {}".format(band_time, args.output, merged))
            shutil.move(merged, args.output)
            os.remove(band_time)
            print("DONE")

        os.remove(tmp_file)