s1base.py

import otbApplication as otb
from Common.otb_numpy_proc import to_otb_pipeline
import os
import glob
import numpy as np
import uuid
import xml.etree.ElementTree as ET
import shutil
from itertools import groupby
from Common.geotools import get_query_bbox
from eodag import EODataAccessGateway

def fetch(shp, dt, output_fld, credentials):
    bbox = get_query_bbox(shp)
    dag = EODataAccessGateway(user_conf_file_path=credentials)
    dag.set_preferred_provider("peps")
    search_criteria = {
        "productType": "S1_SAR_GRD",
        "start": dt.split("/")[0],
        "end": dt.split("/")[1],
        "geom": {"lonmin": bbox[0], "latmin": bbox[1], "lonmax": bbox[2], "latmax": bbox[3]}
    }
    res = dag.search_all(**search_criteria)
    res.filter_property(sensorMode='IW')
    if len(res) > 0:
        os.makedirs(output_fld, exist_ok=True)
        dag.download_all(res, outputs_prefix=output_fld, extract=True)

    # return S1GRDPipeline(output_fld)
    return


class S1GRDPipeline:
    # --- BEGIN SENSOR PROTOTYPE ---
    NAME = 'S1-IW-GRD'
    VAL_TYPE = otb.ImagePixelType_int16
    TMP_TYPE = otb.ImagePixelType_float
    PTRN_dir = 'S1*_IW_GRD*/S1*_IW_GRD*.SAFE'
    PTRN_ref = '-iw-grd-'
    VH_name = 'vh'
    PTRN = ['measurement/s1*-iw-grd-vh-*.tiff', 'measurement/s1*-iw-grd-vv-*.tiff']
    FEAT_exp = {
        'VH': 'im1b1',
        'VV': 'im2b1',
        'VH_db': '1000*log10(abs(im1b1)+1e-6)',
        'VV_db': '1000*log10(abs(im2b1)+1e-6)',
        'POL_RATIO': 'im1b1/im2b1'
    }
    NDT = 0.0

    @classmethod
    def _check(cls, x):
        lst = os.path.basename(cls.PTRN_dir).split('*')
        return all([t in os.path.basename(x) for t in lst])

    @classmethod
    def _img_id(cls, x):
        if cls._check(x):
            return '_'.join(os.path.basename(x).split('_')[4:7])
        else:
            return None

    @classmethod
    def _img_date(cls, x):
        if cls._check(x):
            return os.path.basename(x).split('_')[4].split('T')[0]
        else:
            return None

    @classmethod
    def _check_roi(cls, x, roi, min_surf=0.0):
        if cls._check(x):
            if os.path.isdir(x):
                er = otb.Registry.CreateApplication('ExtractROI')
                bnd = glob.glob(os.path.join(x, cls.PTRN[0]))[0]
                er.SetParameterString('in', bnd)
                er.SetParameterString('mode', 'fit')
                er.SetParameterString('mode.fit.im', roi)
                try:
                    er.Execute()
                except:
                    return False
                arr = er.GetImageAsNumpyArray('out')
                if (np.sum(arr != cls.NDT) / (arr.shape[0] * arr.shape[1])) <= min_surf:
                    return False
                else:
                    return True

    @classmethod
    def _check_direction(cls, x):
        root = ET.parse(os.path.join(x, 'manifest.safe')).getroot()
        for item in root.iter():
            if item.tag.endswith('pass'):
                if item.text == 'ASCENDING':
                    return 0
                elif item.text == 'DESCENDING':
                    return 1
                else:
                    return None

    @classmethod
    def _get_stiched_filename(cls, fn):
        tm1 = os.path.basename(fn).split('-')[-4].split('t')[1]
        tm2 = os.path.basename(fn).split('-')[-5].split('t')[1]
        fn = fn.replace(tm1, 'xxxxxx').replace(tm2, 'xxxxxx')
        return fn

    def _check_satellite(self, x):
        return os.path.basename(x).split('_')[0].lower()

    def __init__(self, fld, roi, temp_fld='/tmp', input_date_interval=None, roi_min_surf=0.0,
                 direction=None, satellite=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.temp_fld = temp_fld + os.sep + self.id
        self.direction = direction
        self.satellite = satellite
        self.roi = roi
        self.image_list = self.parse_folder(self.folder, roi_min_surf)
        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:
            for p in self.PTRN:
                ifn, ofn = self.get_file(img, p)
                self.append(to_otb_pipeline(ifn), ofn, self.VAL_TYPE, 'out', is_output=True)

    def __del__(self):
        if os.path.exists(self.temp_fld):
            shutil.rmtree(self.temp_fld)

    def parse_folder(self, fld, roi_min_surf=0.0):
        img_list = [os.path.abspath(x) for x in glob.glob(os.path.join(fld, self.PTRN_dir))
                    if os.path.isdir(x) and self._check_roi(x, self.roi, roi_min_surf)]

        if self.satellite is not None:
            img_list = [x for x in img_list if self._check_satellite(x) == self.satellite]

        asc, desc = [], []
        for x in img_list:
            if self._check_direction(x) == 0:
                asc.append(x)
            elif self._check_direction(x) == 1:
                desc.append(x)

        out = None
        if self.direction == 'ascending':
            out = asc
        elif self.direction == 'descending':
            out = desc
        elif self.direction == None:
            if len(asc) >= len(desc):
                self.direction = 'ascending'
                out = asc
            else:
                self.direction = 'descending'
                out = desc
            print('[INFO] No direction selected, returning majority : ' + self.direction)
        return sorted(out, key=lambda x: self._img_id(x))

    def get_file(self, img, ptrn):
        in_fn = glob.glob(os.path.join(img, ptrn))[0]
        out_fn = in_fn.replace(self.folder, '').lstrip(os.sep)
        return in_fn, os.path.splitext(out_fn)[0] + '.tiff'

    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 calibrate(self, lut='sigma'):
        proc_idx = self.out_idx.copy()
        self.out_idx = []
        for t in proc_idx:
            sarcal = otb.Registry.CreateApplication('SARCalibration')
            sarcal.SetParameterInputImage('in', self.pipe[t].GetParameterOutputImage(self.out_p[t]))
            sarcal.SetParameterString('lut', lut)
            sarcal.Execute()
            fn = self.files[t].replace('.tiff', '_calib.tiff')
            ty = self.TMP_TYPE
            self.append(sarcal, fn, ty, 'out', is_output=True)

    def orthorectify(self, dem_fld, geoid=None, grid_spacing=40):
        assert(os.path.isdir(dem_fld))
        proc_idx = self.out_idx.copy()
        self.out_idx = []
        for t in proc_idx:
            ortho = otb.Registry.CreateApplication('OrthoRectification')
            ortho.SetParameterString('elev.dem', dem_fld)
            if geoid is not None:
                assert (os.path.exists(geoid))
                ortho.SetParameterString('elev.geoid', geoid)
            ortho.SetParameterInputImage('io.in', self.pipe[t].GetParameterOutputImage(self.out_p[t]))
            ortho.SetParameterInt('opt.gridspacing', grid_spacing)
            ortho.SetParameterString('outputs.mode', 'orthofit')
            ortho.SetParameterString('outputs.ortho', self.roi)
            ortho.Execute()
            fn = self.files[t].replace('.tiff', '_ortho.tiff')
            ty = self.TMP_TYPE
            self.append(ortho, fn, ty, 'io.out', is_output=True)

    def superimpose(self):
        proc_idx = self.out_idx.copy()
        self.out_idx = []
        for t in proc_idx:
            si = otb.Registry.CreateApplication('Superimpose')
            si.SetParameterString('inr', self.roi)
            si.SetParameterInputImage('inm', self.pipe[t].GetParameterOutputImage(self.out_p[t]))
            si.Execute()
            fn = self.files[t].replace('.tiff', '_roi.tiff')
            self.append(si, fn, self.TMP_TYPE, 'out', is_output=True)

    def stitch(self):
        proc_idx = self.out_idx.copy()
        self.out_idx = []
        q = [[self.image_list.index(v) for v in list(i)]
             for j, i in groupby(self.image_list, lambda x: self._img_date(x))]
        for s in q:
            if len(s) == 1:
                self.append(self.pipe[proc_idx[2 * s[0]]], self.files[proc_idx[2 * s[0]]],
                            self.types[proc_idx[2 * s[0]]], self.out_p[proc_idx[2 * s[0]]], is_output=True)
                self.append(self.pipe[proc_idx[2 * s[0] + 1]], self.files[proc_idx[2 * s[0] + 1]],
                            self.types[proc_idx[2 * s[0] + 1]], self.out_p[proc_idx[2 * s[0] + 1]], is_output=True)
            else:
                for k in range(2):
                    fn = self._get_stitched_filename(self.files[proc_idx[2*s[0]+k]])
                    bm = otb.Registry.CreateApplication('BandMathX')
                    [bm.AddImageToParameterInputImageList('il', self.pipe[proc_idx[2*u+k]].GetParameterOutputImage(self.out_p[proc_idx[2*u+k]])) for u in s]
                    bm.SetParameterString('exp',
                                          'vmax({' + ';'.join(['im%db1' % (i + 1) for i in range(len(s))]) + '})')
                    bm.Execute()
                    self.append(bm, fn, self.types[proc_idx[2*s[0]+k]], 'out', is_output=True)

    def multitemp_speckle_filter(self, win_size=3, outcore_on_disk=True):
        proc_idx = self.out_idx.copy()
        self.out_idx = []
        tmp_pipe = []
        tmp_fns = []
        ty = self.TMP_TYPE
        for u in range(2):
            oc = otb.Registry.CreateApplication('MultitempFilteringOutcore')
            for t in proc_idx[u::2]:
                oc.AddImageToParameterInputImageList('inl', self.pipe[t].GetParameterOutputImage(self.out_p[t]))
            oc.SetParameterString('wr', str(win_size))
            if not outcore_on_disk:
                oc.Execute()
                oc_idx = len(self.pipe)
                self.append(oc)
            else:
                if not os.path.exists(self.temp_fld):
                    os.makedirs(self.temp_fld)
                oc_fn = os.path.join(self.temp_fld, 'outcore_{}.tif'.format(['vh', 'vv'][u]))
                oc.SetParameterString('oc', oc_fn)
                oc.ExecuteAndWriteOutput()
            for t in proc_idx[u::2]:
                smooth = otb.Registry.CreateApplication('Smoothing')
                smooth.SetParameterInputImage('in', self.pipe[t].GetParameterOutputImage(self.out_p[t]))
                smooth.SetParameterString('type', 'mean')
                smooth.SetParameterString('type.mean.radius', str(win_size))
                smooth.Execute()
                self.append(smooth)
                bm = otb.Registry.CreateApplication('BandMath')
                if not outcore_on_disk:
                    bm.AddImageToParameterInputImageList('il', self.pipe[oc_idx].GetParameterOutputImage('oc'))
                else:
                    bm.SetParameterStringList('il', [oc_fn])
                bm.AddImageToParameterInputImageList('il', self.pipe[-1].GetParameterOutputImage('out'))
                bm.SetParameterString('exp', 'im2b1*im1b1/im1b2')
                bm.Execute()
                tmp_fns.append(self.files[t].replace('.tiff', '_filt.tiff'))
                tmp_pipe.append(bm)

        N = int(len(tmp_pipe) / 2)
        for i in range(N):
            self.append(tmp_pipe[i], tmp_fns[i], ty, 'out', is_output=True)
            self.append(tmp_pipe[N + i], tmp_fns[N + i], ty, 'out', is_output=True)

    def compute_features(self, feat_list=['VH_db', 'VV_db']):
        proc_idx = self.out_idx.copy()
        self.out_idx = []
        for k in range(0,len(proc_idx),2):
            bm = otb.Registry.CreateApplication('BandMathX')
            bm.AddImageToParameterInputImageList('il', self.pipe[proc_idx[k]].GetParameterOutputImage(
                self.out_p[proc_idx[k]]))
            bm.AddImageToParameterInputImageList('il', self.pipe[proc_idx[k+1]].GetParameterOutputImage(
                self.out_p[proc_idx[k+1]]))
            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('-vh', '-feat')
            if all(['db' in x for x in feat_list]):
                ty = self.VAL_TYPE
            else:
                ty = self.TMP_TYPE
            self.append(bm, fn, ty, '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:
            er = otb.Registry.CreateApplication('ExtractROI')
            er.SetParameterInputImage('in', self.pipe[t].GetParameterOutputImage(self.out_p[t]))
            er.SetParameterString('mode', 'fit')
            er.SetParameterString('mode.fit.im', self.roi)
            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))
            er.SetParameterString('out', out_file)
            er.SetParameterOutputImagePixelType('out', self.types[t])
            er.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