planet_mosaics.py

import otbApplication as otb
from Common.otb_numpy_proc import to_otb_pipeline
from osgeo import ogr
from pyproj import Transformer as T
import requests
import os
from tqdm import tqdm
import urllib.request
from datetime import datetime
from dateutil.relativedelta import *
import uuid
import glob

def fetch(shp, dt, out_fld, api_key):
    ds = ogr.Open(shp)
    ly = ds.GetLayer(0)
    xt = ly.GetExtent()
    shp_srs = int(ly.GetSpatialRef().GetAuthorityCode(None))
    qry_srs = 4326

    i_bbox = [xt[0], xt[2], xt[1], xt[3]]
    if shp_srs == qry_srs:
        bbox = i_bbox
    else:
        tr = T.from_crs(shp_srs, qry_srs, always_xy=True)
        lon_min, lat_min = tr.transform(xt[0], xt[2])
        lon_max, lat_max = tr.transform(xt[1], xt[3])
        bbox = [lon_min, lat_min, lon_max, lat_max]

    str_bbox = ','.join(map(str, bbox))

    API_URL = "https://api.planet.com/basemaps/v1/mosaics"
    session = requests.Session()
    session.auth = (api_key, "")

    # date format "YYYY-MM/YYYY-MM"
    sd = datetime.strptime(dt.split('/')[0], '%Y-%m')
    ed = datetime.strptime(dt.split('/')[1], '%Y-%m')
    delta = relativedelta(ed, sd)
    mosaic_dates = [str(d.year) + '-' + str(d.month).zfill(2)
                    for d in (sd + relativedelta(months=+n) for n in range(delta.months+1))]

    if not os.path.exists(out_fld):
        os.makedirs(out_fld)

    for md in mosaic_dates:
        parameters = {
            "name__is": "planet_medres_normalized_analytic_{}_mosaic".format(md)
        }
        res = session.get(API_URL, params=parameters)
        mosaic = res.json()
        if len(mosaic['mosaics']) > 0:
            mosaic_id = mosaic['mosaics'][0]['id']

            search_parameters = {
                'bbox': str_bbox,
                'minimal': True
            }
            quads_url = "{}/{}/quads".format(API_URL, mosaic_id)
            res = session.get(quads_url, params=search_parameters, stream=True)
            quads = res.json()

            l_quads = []
            for q in tqdm(quads['items'], desc='Downloading quads for {}'.format(md)):
                link = q['_links']['download']
                l_quads.append(os.path.join(out_fld, q['id'] + '.tif'))
                if not os.path.exists(l_quads[-1]):
                    urllib.request.urlretrieve(link, l_quads[-1])

            app_mos = otb.Registry.CreateApplication('Mosaic')
            app_mos.SetParameterStringList('il', l_quads)
            app_mos.Execute()
            app_roi = otb.Registry.CreateApplication('ExtractROI')
            app_roi.SetParameterInputImage('in', app_mos.GetParameterOutputImage('out'))
            fn = 'PlanetMosaic_{}_{}.tif'.format(os.path.splitext(os.path.basename(shp))[0], md)
            app_roi.SetParameterString('mode', 'fit')
            app_roi.SetParameterString('mode.fit.vect', shp)
            app_roi.SetParameterString('out', os.path.join(out_fld, fn))
            app_roi.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint16)
            app_roi.ExecuteAndWriteOutput()

            [os.remove(f) for f in l_quads]

    return PlanetMosaicPipeline(out_fld)

class PlanetMosaicPipeline:
    # --- BEGIN SENSOR PROTOTYPE ---

    NAME = 'PLANET-MOSAICS'
    REF_TYPE = otb.ImagePixelType_uint16
    FEAT_TYPE = otb.ImagePixelType_float
    PTRN_name = 'PlanetMosaic_*.tif'
    FEAT_exp = {
        'B': 'im1b1',
        'G': 'im1b2',
        'R': 'im1b3',
        'NIR': 'im1b4',
        'NDVI': '(im1b4-im1b3)/(im1b4+im1b3+1e-6)'
    }
    NDT = 0

    # ---- END SENSOR PROTOTYPE ----

    @classmethod
    def _check(cls,x):
        return (cls.PTRN_name.split('*')[0] in os.path.basename(x) and
                cls.PTRN_name.split('*')[1] in os.path.basename(x))

    @classmethod
    def _img_id(cls,x):
        if cls._check(x):
            return os.path.splitext(os.path.basename(x))
        else:
            return None

    @classmethod
    def _img_date(cls,x):
        if cls._check(x):
            return os.path.splitext(os.path.basename(x))[0].split('_')[-1]
        else:
            return None
    def __init__(self, fld, input_date_interval=None):
        self.pipe = []
        self.files = []
        self.types = []
        self.out_p = []
        self.out_idx = []

        self.id = str(uuid.uuid4())
        self.folder = os.path.abspath(fld)
        self.image_list = self.parse_folder(self.folder)
        self.input_dates = [self._img_date(x) for x in self.image_list]

        if input_date_interval is not None:
            idx = [i for i in range(len(self.input_dates)) if
                   input_date_interval[0] <= self.input_dates[i] <= input_date_interval[1]]
            self.image_list = [self.image_list[i] for i in idx]
            self.input_dates = [self.input_dates[i] for i in idx]

        for img in self.image_list:
            self.append(to_otb_pipeline(img), os.path.basename(img), self.REF_TYPE, 'out', is_output=True)

    def __del__(self):
        return

    def parse_folder(self, fld):
        img_list = [os.path.abspath(x) for x in glob.glob(os.path.join(fld, self.PTRN_name))]
        return sorted(img_list, key=lambda x: self._img_id(x))

    def append(self, app, fname=None, ftype=None, outp=None, is_output=False):
        if is_output:
            self.out_idx.append(len(self.pipe))
        self.pipe.append(app)
        self.files.append(fname)
        self.types.append(ftype)
        self.out_p.append(outp)

    def compute_features(self, feat_list=['B', 'G', 'R', 'NIR', 'NDVI']):
        proc_idx = self.out_idx.copy()
        self.out_idx = []
        for k in range(len(proc_idx)):
            bm = otb.Registry.CreateApplication('BandMathX')
            bm.AddImageToParameterInputImageList('il', self.pipe[k].GetParameterOutputImage(self.out_p[k]))
            expr = []
            for f in feat_list:
                expr.append(self.FEAT_exp[f])
            expr = '{' + ';'.join(expr) + '}'
            bm.SetParameterString('exp', expr)
            bm.Execute()
            fn = self.files[proc_idx[k]].replace('.tif', '_FEAT.tif')
            self.append(bm, fn, self.FEAT_TYPE, 'out', is_output=True)

    def write_outputs(self, fld, update_pipe=False, compress=False):
        out = []
        proc_idx = self.out_idx.copy()
        if update_pipe:
            self.out_idx = []
        for t in proc_idx:
            out_file = os.path.join(fld, self.files[t])
            if compress:
                out_file += '?gdal:co:compress=deflate'
            if not os.path.exists(os.path.dirname(out_file)):
                os.makedirs(os.path.dirname(out_file))
            self.pipe[t].SetParameterString(self.out_p[t], out_file)
            self.pipe[t].SetParameterOutputImagePixelType(self.out_p[t], self.types[t])
            self.pipe[t].ExecuteAndWriteOutput()
            out.append(out_file)
            if update_pipe:
                self.append(to_otb_pipeline(out_file), self.files[t], self.types[t], 'out', is_output=True)
        return out