diff --git a/Common/otb_numpy_proc.py b/Common/otb_numpy_proc.py
new file mode 100644
index 0000000000000000000000000000000000000000..1465d0402ee9d60f825e514df165cb2a64908099
--- /dev/null
+++ b/Common/otb_numpy_proc.py
@@ -0,0 +1,221 @@
+import otbApplication as otb
+import numpy as np
+import sys
+from time import sleep
+import os, psutil
+from functools import partial
+import string, random
+from gdal import BuildVRT
+
+def randomword(length):
+ return ''.join(random.choice(string.ascii_lowercase) for i in range(length))
+
+def generate_stream_rois(img_size, n_lines):
+ streams = []
+ n_streams = int(img_size[1]/n_lines)
+ for i in range(n_streams):
+ stream = {}
+ stream['sizex'] = img_size[0]
+ stream['sizey'] = n_lines
+ stream['startx'] = 0
+ stream['starty'] = i * n_lines
+ streams.append(stream)
+
+ if (img_size[1] % n_lines) :
+ stream = {}
+ stream['sizex'] = img_size[0]
+ stream['sizey'] = img_size[1] % n_lines
+ stream['startx'] = 0
+ stream['starty'] = n_streams * n_lines
+ streams.append(stream)
+
+ return streams
+
+def pad_format_tile(tile, img_size, pad):
+ new_tile = np.array(tile) + np.array([-pad, -pad, pad, pad])
+ pad_out = [0,0,0,0]
+
+ if new_tile[0] < 0:
+ pad_out[0] = pad
+ new_tile[0] = 0
+ if new_tile[1] < 0:
+ pad_out[1] = pad
+ new_tile[1] = 0
+ if new_tile[2] > img_size[0]-1:
+ pad_out[2] = pad
+ new_tile[2] = img_size[0]-1
+ if new_tile[3] > img_size[1]-1:
+ pad_out[3] = pad
+ new_tile[3] = img_size[1]-1
+
+ new_tile[2] = new_tile[2] - new_tile[0] + 1
+ new_tile[3] = new_tile[3] - new_tile[1] + 1
+ pads1 = (pad_out[0], pad_out[2])
+ pads2 = (pad_out[1], pad_out[3])
+
+ return ([int(x) for x in new_tile],pads1,pads2)
+
+
+def generate_tile_rois(img_size, tile_size, pad = 0):
+ tiles = []
+ n_tiles = [int(img_size[0]/tile_size[0]), int(img_size[1]/tile_size[1])]
+ for i in range(2):
+ if img_size[i] % tile_size[i] :
+ n_tiles[i] += 1
+ for c in range(n_tiles[0]):
+ for r in range(n_tiles[1]):
+ pad_out = [0,0,0,0]
+ bbox = [c*tile_size[0],r*tile_size[1],min(c*tile_size[0]+tile_size[0]-1,img_size[0]-1),min(r*tile_size[1]+tile_size[1]-1,img_size[1]-1)]
+ tiles.append(pad_format_tile(bbox,img_size,pad))
+ return tiles
+
+'''
+def generate_stream_regions(img_size, n_lines):
+ streams = []
+ n_streams = int(img_size[1]/n_lines)
+ for i in range(n_streams):
+ stream = otb.itkRegion()
+ stream['size'][0] = img_size[0]
+ stream['size'][1] = n_lines
+ stream['index'][0] = 0
+ stream['index'][1] = i * n_lines
+ streams.append(stream)
+
+ if (img_size[1] % n_lines) :
+ stream = otb.itkRegion()
+ stream['size'][0] = img_size[0]
+ stream['size'][1] = img_size[1] % n_lines
+ stream['index'][0] = 0
+ stream['index'][1] = n_streams * n_lines
+ streams.append(stream)
+
+ return streams
+'''
+
+'''
+def generate_streams(otb_pipeline):
+ otb_pipeline.Execute()
+ img_size = tuple(otb_pipeline.GetImageSize("out"))
+ stream_rois = generate_stream_rois(img_size,1000)
+ ext_roi_apps = []
+ for roi in stream_rois :
+ print("setting roi : " + str(roi))
+ app = otb.Registry.CreateApplication('ExtractROI')
+ app.SetParameterInputImage("in", otb_pipeline.GetParameterOutputImage("out"))
+ app.SetParameterInt("sizex", roi['sizex'])
+ app.SetParameterInt("sizey", roi['sizey'])
+ app.SetParameterInt("startx", roi['startx'])
+ app.SetParameterInt("starty", roi['starty'])
+ ext_roi_apps.append(app)
+ ext_roi_apps[-1].Execute()
+
+ return ext_roi_apps
+'''
+
+def stream_function(otb_in, function, tile_size, write_output_to=None, pad=None, out_key = 'out', work_dir='.', prefix = None):
+ otb_pipeline = to_otb_pipeline(otb_in)
+
+ if prefix is None:
+ prefix = randomword(16)
+
+ if pad is None:
+ pad = check_function_padding(function)
+
+ otb_pipeline.Execute()
+ img_size = tuple(otb_pipeline.GetImageSize(out_key))
+
+ streams = generate_tile_rois(img_size,tile_size,pad)
+
+ t = 1
+ output_streams = []
+ for stream in streams:
+ roi = otb.Registry.CreateApplication('ExtractROI')
+ roi.SetParameterInputImage('in', otb_pipeline.GetParameterOutputImage(out_key))
+ roi.SetParameterInt('startx', stream[0][0])
+ roi.SetParameterInt('starty', stream[0][1])
+ roi.SetParameterInt('sizex', stream[0][2])
+ roi.SetParameterInt('sizey', stream[0][3])
+ roi.Execute()
+ img = roi.ExportImage('out')
+
+ # Launch numpy function !
+ padw = (stream[1],stream[2]) if img['array'].ndim == 2 else (stream[1],stream[2],(0,0))
+ img['array'] = np.ascontiguousarray(function(np.pad(img['array'],padw)))
+
+ cpy = otb.Registry.CreateApplication('ExtractROI')
+ cpy.ImportVectorImage('in', img)
+ output_streams.append(os.path.join(work_dir, prefix + '_' + str(t) + '.tif'))
+ cpy.SetParameterString('out', output_streams[-1])
+ cpy.ExecuteAndWriteOutput()
+ print(cpy.GetImageSize('out'))
+ app = None
+ img = None
+ t += 1
+
+ vrt = os.path.join(work_dir, prefix + '.vrt')
+ BuildVRT(vrt, output_streams)
+ output_streams.append(vrt)
+
+ new_pipeline = otb.Registry.CreateApplication('ExtractROI')
+ new_pipeline.SetParameterString('in', vrt)
+
+ if write_output_to is not None:
+ new_pipeline.SetParameterString('out', write_output_to)
+ new_pipeline.ExecuteAndWriteOutput()
+ [os.remove(f) for f in output_streams]
+ return write_output_to, []
+ else:
+ return new_pipeline, output_streams
+
+'''
+def merge_streams(ext_roi_apps):
+ app = otb.Registry.CreateApplication('Mosaic')
+ i=0
+ for roi in ext_roi_apps:
+ img = roi.ExportImage('out')
+ img['array'] /= 10
+ print(psutil.Process(os.getpid()).memory_info().rss)
+ sleep(5)
+ app.ImportVectorImage('il',img,i)
+ i += 1
+ return app
+'''
+
+def test_function(arr):
+ #return arr[2:-2,2:-2,:]
+ return np.sum(arr,axis=2)[6:-6,6:-6]
+
+def check_function_padding(function, test_tile_size=(100,100)):
+ try:
+ out = function(np.ones((test_tile_size)))
+ except:
+ out = function(np.ones((test_tile_size[0],test_tile_size[1],3)))
+ return int(max(abs(test_tile_size[0] - out.shape[0]), abs(test_tile_size[1] - out.shape[1])) / 2)
+
+def to_otb_pipeline(obj):
+ if isinstance(obj, str):
+ ppl = otb.Registry.CreateApplication('ExtractROI')
+ ppl.SetParameterString('in',obj)
+ ppl.Execute()
+ elif isinstance(obj,otb.Application):
+ ppl = obj
+ elif isinstance(obj,list): # to-do, check if all sublists are lists of otbApp
+ ppl = obj
+ else:
+ sys.exit("Impossible to convert object ot OTB pipeline")
+ return ppl
+
+def do_something(fn):
+ # Create a smoothing application
+ app = otb.Registry.CreateApplication("ExtractROI")
+ app.SetParameterString("in", fn)
+
+ app2, to_del = stream_function(app, partial(test_function), tile_size=(1000,1000))
+
+ app2.SetParameterString('out','pippo.tif')
+ app2.ExecuteAndWriteOutput()
+
+ [os.remove(f) for f in to_del]
+
+if __name__ == "__main__":
+ do_something("/DATA/Koumbia/Tong/SENTINEL2_Aug_Oct_2017_RGB_NIR_NDVI.tif")
diff --git a/Common/otb_vector_proc.py b/Common/otb_vector_proc.py
new file mode 100644
index 0000000000000000000000000000000000000000..4fd5e8c53656277692851e6a5252e179f38dcba8
--- /dev/null
+++ b/Common/otb_vector_proc.py
@@ -0,0 +1,61 @@
+import gdal, ogr, osr
+import numpy as np
+import otb_numpy_proc
+import otbApplication as otb
+import os
+
+def pixels_to_vector(selection_array, otb_pipeline, out_vector, out_key='out', var_dict=None):
+ geo_transform = otb_pipeline.GetImageMetaData(out_key)['GeoTransform']
+ epsg = int(otb_pipeline.GetImageMetaData(out_key)['ProjectionRef'].split('"')[-2])
+ srs = osr.SpatialReference()
+ srs.ImportFromEPSG(epsg)
+
+ drv = ogr.GetDriverByName('ESRI Shapefile')
+ ds = drv.CreateDataSource(out_vector)
+ ly = ds.CreateLayer(os.path.basename(os.path.splitext(out_vector)[0]), geom_type=ogr.wkbPoint, srs=srs)
+ ly.CreateField(ogr.FieldDefn('id', ogr.OFTInteger64))
+ ly.CreateField(ogr.FieldDefn('px', ogr.OFTInteger64))
+ ly.CreateField(ogr.FieldDefn('py', ogr.OFTInteger64))
+ if var_dict is not None:
+ for v in var_dict:
+ if 'float' in str(v['array'].dtype):
+ typ = ogr.OFTReal
+ else:
+ typ = ogr.OFTInteger
+ if v['array'].ndim == 2:
+ for t in range(v['array'].shape[1]):
+ ly.CreateField(ogr.FieldDefn(v['prefix']+'_%d' % t, typ))
+ else:
+ ly.CreateField(ogr.FieldDefn(v['prefix'], typ))
+
+ i = 0
+ pix_pnts = np.where(selection_array)
+ for p in zip(pix_pnts[0]+0.5,pix_pnts[1]+0.5):
+ X = geo_transform[0] + p[1] * geo_transform[1]
+ Y = geo_transform[3] + p[0] * geo_transform[5]
+ gp = ogr.Geometry(ogr.wkbPoint)
+ gp.AddPoint(X,Y)
+ gf = ogr.Feature(ly.GetLayerDefn())
+ gf.SetGeometry(gp)
+ gf.SetField('id', i)
+ gf.SetField('px', p[1])
+ gf.SetField('py', p[0])
+ if var_dict is not None:
+ for v in var_dict:
+ if v['array'].ndim == 2:
+ for t in range(v['array'].shape[1]):
+ gf.SetField(v['prefix']+'_%d' % t, float(v['array'][i,t]))
+ else:
+ gf.SetField(v['prefix'], float(v['array'][i]))
+ ly.CreateFeature(gf)
+ i += 1
+
+ ds = None
+
+if __name__ == '__main__':
+ p = otb_numpy_proc.to_otb_pipeline('./pippo.tif')
+ p.Execute()
+ arr = np.random.random((100, 100)) < 0.1
+ drv = ogr.GetDriverByName('ESRI Shapefile')
+ drv.DeleteDataSource('./pippuzzo.shp')
+ pixels_to_vector(arr, p, 'pippuzzo.shp')
diff --git a/TimeSeries/s1base.py b/TimeSeries/s1base.py
new file mode 100644
index 0000000000000000000000000000000000000000..e8802bf454480a3fc8ae336b9aa5b8c8bc24266b
--- /dev/null
+++ b/TimeSeries/s1base.py
@@ -0,0 +1,307 @@
+import otbApplication as otb
+from 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
+
+class S1GRDPipeline:
+ # --- BEGIN SENSOR PROTOTYPE ---
+ NAME = 'S1-IW-GRD'
+ VAL_TYPE = otb.ImagePixelType_int16
+ TMP_TYPE = otb.ImagePixelType_float
+ PTRN_dir = '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(x.split('_')[4:7])
+ else:
+ return None
+
+ @classmethod
+ def _img_date(cls, x):
+ if cls._check(x):
+ return 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
+
+ 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):
+ print(glob.glob(os.path.join(fld, self.PTRN_dir)))
+ 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 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.files[proc_idx[2*s[0]+k]]
+ 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')
+ 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
+
+
+
+
+
+
diff --git a/TimeSeries/s2sen2cor.py b/TimeSeries/s2sen2cor.py
new file mode 100644
index 0000000000000000000000000000000000000000..01de7aca13e6ac0886626b55b0d32315e97b3a86
--- /dev/null
+++ b/TimeSeries/s2sen2cor.py
@@ -0,0 +1,173 @@
+import gdal
+import otbApplication as otb
+from otb_numpy_proc import to_otb_pipeline
+import glob
+import os
+import sys
+import zipfile
+import osr
+import datetime
+from s2theia import *
+
+class S2Sen2CorTilePipeline(S2TheiaTilePipeline):
+ # --- BEGIN SENSOR PROTOTYPE ---
+
+ NAME = 'S2-SEN2COR'
+ REF_TYPE = otb.ImagePixelType_uint16
+ MSK_TYPE = otb.ImagePixelType_uint8
+ PTRN_dir = 'S2*_MSIL2A*'
+ PTRN_ref = '_B'
+ B2_name = 'B02'
+ PTRN_10m = ['GRANULE/*/IMG_DATA/R10m/*_B02_10m.jp2',
+ 'GRANULE/*/IMG_DATA/R10m/*_B03_10m.jp2',
+ 'GRANULE/*/IMG_DATA/R10m/*_B04_10m.jp2',
+ 'GRANULE/*/IMG_DATA/R10m/*_B08_10m.jp2']
+ PTRN_20m = ['GRANULE/*/IMG_DATA/R20m/*_B05_20m.jp2',
+ 'GRANULE/*/IMG_DATA/R20m/*_B06_20m.jp2',
+ 'GRANULE/*/IMG_DATA/R20m/*_B07_20m.jp2',
+ 'GRANULE/*/IMG_DATA/R20m/*_B8A_20m.jp2',
+ 'GRANULE/*/IMG_DATA/R20m/*_B11_20m.jp2',
+ 'GRANULE/*/IMG_DATA/R20m/*_B12_20m.jp2']
+ PTRN_msk = ['GRANULE/*/IMG_DATA/R20m/*_SCL_20m.jp2']
+ PTRN_ful = PTRN_10m[0:3] + PTRN_20m[0:3] + [PTRN_10m[3]] + PTRN_20m[3:]
+ FEAT_exp = {
+ 'B2': 'im1b1',
+ 'B3': 'im1b2',
+ 'B4': 'im1b3',
+ 'B5': 'im1b4',
+ 'B6': 'im1b5',
+ 'B7': 'im1b6',
+ 'B8': 'im1b7',
+ 'B8A': 'im1b8',
+ 'B11': 'im1b9',
+ 'B12': 'im1b10',
+ 'NDVI': '(im1b7-im1b3)/(im1b7+im1b3+1e-6)',
+ 'NDWI': '(im1b2-im1b7)/(im1b2+im1b7+1e-6)',
+ 'BRI': 'sqrt(' + '+'.join(['im1b%d*im1b%d' % (i, i) for i in range(1, 11)]) + ')',
+ 'MNDWI': '(im1b2-im1b9)/(im1b2+im1b9+1e-6)',
+ 'SWNDVI': '(im1b9-im1b7)/(im1b9+im1b7+1e-6)',
+ 'NDRE': '(im1b7-im1b4)/(im1b7+im1b4+1e-6)'
+ }
+ NDT = 0
+
+ @classmethod
+ def _check(cls,x):
+ return cls.PTRN_dir.split('_')[-1].replace('*', '') in os.path.basename(x)
+
+ @classmethod
+ def _img_id(cls,x):
+ if cls._check(x):
+ if os.path.isdir(x) or os.path.splitext(x)[-1] == '.zip':
+ return '_'.join(os.path.splitext(x)[0].split('_')[-5:])
+ else:
+ return None
+ else:
+ return None
+
+ @classmethod
+ def _img_date(cls,x):
+ if cls._check(x):
+ if os.path.isdir(x) or os.path.splitext(x)[-1] == '.zip':
+ return os.path.splitext(x)[0].split('_')[-5].split('T')[0]
+ else:
+ return None
+ else:
+ return None
+
+ @classmethod
+ def _tile_id(cls,x):
+ if cls._check(x):
+ if os.path.isdir(x) or os.path.splitext(x)[-1] == '.zip':
+ return os.path.splitext(x)[0].split('_')[-2]
+ else:
+ return None
+ else:
+ return None
+
+ @classmethod
+ def _tile_cloud_percentage(cls, x):
+ if cls._check(x):
+ if os.path.isdir(x):
+ fid = open(glob.glob(os.path.join(x, '*/*/MTD_TL.xml'))[0], 'r')
+ mtd = fid.read()
+ elif os.path.splitext(x)[-1] == '.zip':
+ arch = zipfile.ZipFile(x)
+ fid = [name for name in arch.namelist() if name.endswith('MTD_TL.xml')]
+ mtd = arch.read(fid[0])
+ root = et.fromstring(mtd)
+ f = filter(lambda x: x.tag == "CLOUDY_PIXEL_PERCENTAGE", root.findall('*/*/*'))
+ r = list(f)
+ return float(r[0].text)
+
+ """
+ def _process_mask(self, msks):
+ msk_pipe = [otb.Registry.CreateApplication('Superimpose')]
+ msk_pipe[-1].SetParameterInputImage('inr', self.pipe[-1].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterInputImage('inm', msks[0].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterString('interpolator', 'nn')
+ msk_pipe[-1].Execute()
+ msk_pipe.append(otb.Registry.CreateApplication('BandMath'))
+ msk_pipe[-1].AddImageToParameterInputImageList('il', msk_pipe[-2].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterString('exp', 'im1b1==0 || im1b1==1 || im1b1==2 || im1b1==7 ||im1b1==3 || im1b1==8 || im1b1==9 || im1b1==10')
+ msk_pipe[-1].Execute()
+ msk_pipe.append(otb.Registry.CreateApplication('BinaryMorphologicalOperation'))
+ msk_pipe[-1].SetParameterInputImage('in', msk_pipe[-2].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterString('filter', 'erode')
+ msk_pipe[-1].SetParameterString('structype', 'ball')
+ msk_pipe[-1].SetParameterString('xradius', '5')
+ msk_pipe[-1].SetParameterString('yradius', '5')
+ msk_pipe[-1].Execute()
+ msk_pipe.append(otb.Registry.CreateApplication('BinaryMorphologicalOperation'))
+ msk_pipe[-1].SetParameterInputImage('in', msk_pipe[-2].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterString('filter', 'dilate')
+ msk_pipe[-1].SetParameterString('structype', 'ball')
+ msk_pipe[-1].SetParameterString('xradius', '20')
+ msk_pipe[-1].SetParameterString('yradius', '20')
+ msk_pipe[-1].Execute()
+ return msk_pipe
+ """
+
+ def _process_mask(self, msks):
+ msk_pipe = [otb.Registry.CreateApplication('BandMath')]
+ msk_pipe[-1].AddImageToParameterInputImageList('il', msks[0].GetParameterOutputImage('out'))
+ # excluding thin cirrus from masking : no im1b1 == 10 below
+ # excluding cloud shadows here (for opening) : no im1b1 == 3 below
+ msk_pipe[-1].SetParameterString('exp', 'im1b1==0 || im1b1==1 || im1b1==8 || im1b1==9')
+ msk_pipe[-1].Execute()
+ msk_pipe.append(otb.Registry.CreateApplication('BandMath'))
+ msk_pipe[-1].AddImageToParameterInputImageList('il', msks[0].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterString('exp', 'im1b1==2 || im1b1 == 3 || im1b1==7')
+ msk_pipe[-1].Execute()
+ msk_pipe.append(otb.Registry.CreateApplication('BinaryMorphologicalOperation'))
+ msk_pipe[-1].SetParameterInputImage('in', msk_pipe[-2].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterString('filter', 'opening')
+ msk_pipe[-1].SetParameterString('structype', 'ball')
+ msk_pipe[-1].SetParameterString('xradius', '5')
+ msk_pipe[-1].SetParameterString('yradius', '5')
+ msk_pipe[-1].Execute()
+ msk_pipe.append(otb.Registry.CreateApplication('BandMath'))
+ msk_pipe[-1].AddImageToParameterInputImageList('il', msk_pipe[-4].GetParameterOutputImage('out'))
+ msk_pipe[-1].AddImageToParameterInputImageList('il', msk_pipe[-2].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterString('exp', 'im1b1 || im2b1')
+ msk_pipe[-1].Execute()
+ msk_pipe.append(otb.Registry.CreateApplication('BinaryMorphologicalOperation'))
+ msk_pipe[-1].SetParameterInputImage('in', msk_pipe[-2].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterString('filter', 'dilate')
+ msk_pipe[-1].SetParameterString('structype', 'ball')
+ msk_pipe[-1].SetParameterString('xradius', '20')
+ msk_pipe[-1].SetParameterString('yradius', '20')
+ msk_pipe[-1].Execute()
+ msk_pipe.append(otb.Registry.CreateApplication('Superimpose'))
+ msk_pipe[-1].SetParameterInputImage('inr', self.pipe[-1].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterInputImage('inm', msk_pipe[-2].GetParameterOutputImage('out'))
+ msk_pipe[-1].SetParameterString('interpolator', 'nn')
+ msk_pipe[-1].Execute()
+ return msk_pipe
+
+ # ---- END SENSOR PROTOTYPE ----
+
+class S2Sen2CorPipeline(S2TheiaPipeline):
+
+ S2TilePipeline = S2Sen2CorTilePipeline
+ _check = S2TilePipeline._check
+ _tile_id = S2TilePipeline._tile_id
diff --git a/TimeSeries/s2theia.py b/TimeSeries/s2theia.py
new file mode 100644
index 0000000000000000000000000000000000000000..9d2b7d940b63e89d237ecca5bb4074fc101f478c
--- /dev/null
+++ b/TimeSeries/s2theia.py
@@ -0,0 +1,810 @@
+import sys
+import warnings
+
+import gdal
+import otbApplication as otb
+
+from otb_numpy_proc import to_otb_pipeline
+import numpy as np
+import glob
+import os
+import xml.etree.ElementTree as et
+import zipfile
+import osr
+import datetime
+import uuid
+import shutil
+
+class S2TheiaTilePipeline:
+ # --- BEGIN SENSOR PROTOTYPE ---
+
+ NAME = 'S2-THEIA'
+ REF_TYPE = otb.ImagePixelType_int16
+ MSK_TYPE = otb.ImagePixelType_uint8
+ PTRN_dir = 'SENTINEL2*'
+ PTRN_ref = '_FRE_'
+ B2_name = 'B2'
+ PTRN_10m = ['*_FRE_B2.tif', '*_FRE_B3.tif', '*_FRE_B4.tif', '*_FRE_B8.tif']
+ PTRN_20m = ['*_FRE_B5.tif', '*_FRE_B6.tif', '*_FRE_B7.tif', '*_FRE_B8A.tif', '*_FRE_B11.tif', '*_FRE_B12.tif']
+ PTRN_msk = ['MASKS/*_EDG_R1.tif', 'MASKS/*_SAT_R1.tif', 'MASKS/*_CLM_R1.tif']
+ PTRN_ful = PTRN_10m[0:3] + PTRN_20m[0:3] + [PTRN_10m[3]] + PTRN_20m[3:]
+ FEAT_exp = {
+ 'B2': 'im1b1',
+ 'B3': 'im1b2',
+ 'B4': 'im1b3',
+ 'B5': 'im1b4',
+ 'B6': 'im1b5',
+ 'B7': 'im1b6',
+ 'B8': 'im1b7',
+ 'B8A': 'im1b8',
+ 'B11': 'im1b9',
+ 'B12': 'im1b10',
+ 'NDVI': '(im1b7-im1b3)/(im1b7+im1b3+1e-6)',
+ 'NDWI': '(im1b2-im1b7)/(im1b2+im1b7+1e-6)',
+ 'BRI': 'sqrt(' + '+'.join(['im1b%d*im1b%d' % (i, i) for i in range(1, 11)]) + ')',
+ 'MNDWI': '(im1b2-im1b9)/(im1b2+im1b9+1e-6)',
+ 'SWNDVI': '(im1b9-im1b7)/(im1b9+im1b7+1e-6)',
+ 'NDRE': '(im1b7-im1b4)/(im1b7+im1b4+1e-6)'
+ }
+ NDT = -10000
+
+ @classmethod
+ def _check(cls,x):
+ return cls.PTRN_dir.replace('*', '') in os.path.basename(x)
+
+ @classmethod
+ def _img_id(cls,x):
+ if cls._check(x):
+ if os.path.isdir(x):
+ return x.split('_')[-5]
+ elif os.path.splitext(x)[-1] == '.zip':
+ return x.split('_')[-4]
+ else:
+ return None
+ else:
+ return None
+
+ @classmethod
+ def _img_date(cls,x):
+ if cls._check(x):
+ if os.path.isdir(x):
+ return x.split('_')[-5].split('-')[0]
+ elif os.path.splitext(x)[-1] == '.zip':
+ return x.split('_')[-4].split('-')[0]
+ else:
+ return None
+ else:
+ return None
+
+ @classmethod
+ def _tile_id(cls,x):
+ if cls._check(x):
+ if os.path.isdir(x):
+ return x.split('_')[-3]
+ elif os.path.splitext(x)[-1] == '.zip':
+ return x.split('_')[-2]
+ else:
+ return None
+ else:
+ return None
+
+ @classmethod
+ def _tile_cloud_percentage(cls, x):
+ if cls._check(x):
+ if os.path.isdir(x):
+ fid = open(glob.glob(os.path.join(x, '*_MTD_ALL.xml'))[0], 'r')
+ mtd = fid.read()
+ elif os.path.splitext(x)[-1] == '.zip':
+ arch = zipfile.ZipFile(x)
+ fid = [name for name in arch.namelist() if name.endswith('_MTD_ALL.xml')]
+ mtd = arch.read(fid[0])
+ root = et.fromstring(mtd)
+ f = filter(lambda x: x.get('name') == 'CloudPercent', root.findall('*/*/*/*/*/*'))
+ r = list(f)
+ return float(r[0].text)
+
+ @classmethod
+ def _check_roi(cls, x, roi, min_surf=0.0, temp_fld='/tmp'):
+ if cls._check(x):
+ if os.path.isdir(x):
+ er = otb.Registry.CreateApplication('ExtractROI')
+ bnd = glob.glob(os.path.join(x, cls.PTRN_20m[0]))[0]
+ elif os.path.splitext(x)[-1] == '.zip':
+ idf = cls.PTRN_20m[0].replace('*', '')
+ arch = zipfile.ZipFile(x)
+ fid = [name for name in arch.namelist() if idf in name]
+ fle = arch.read(fid[0])
+ bnd = os.path.join(temp_fld, os.path.basename(fid[0]))
+ tgt = open(bnd, 'wb')
+ with fle,tgt:
+ shutil.copyfileobj(fle, tgt)
+ er.SetParameterString('in', bnd)
+ er.SetParameterString('mode', 'fit')
+ er.SetParameterString('mode.fit.vect', roi)
+ er.Execute()
+ arr = er.GetImageAsNumpyArray('out')
+ if (np.sum(arr != cls.NDT) / (arr.shape[0]*arr.shape[1])) <= min_surf:
+ return False
+ else:
+ return True
+
+ def _process_mask(self, msks):
+ msk_pipe = [otb.Registry.CreateApplication('BandMath')]
+ [msk_pipe[-1].AddImageToParameterInputImageList('il', x.GetParameterOutputImage('out')) for x in
+ msks]
+ msk_pipe[-1].SetParameterString('exp', 'im1b1!=0 || im2b1!=0 || im3b1!=0')
+ msk_pipe[-1].Execute()
+ return msk_pipe
+
+ # ---- END SENSOR PROTOTYPE ----
+
+ def __init__(self, fld, tile, temp_fld='/tmp', input_date_interval=None, max_clouds_percentage=None,
+ filter_by_roi=None, roi_min_surf=0.0, dummy_read=False):
+ self.pipe = []
+ self.files = []
+ self.types = []
+ self.out_p = []
+ self.out_idx = []
+
+ self.id = str(uuid.uuid4())
+ self.folder = os.path.abspath(fld)
+ self.tile_id = tile
+ self.temp_fld = temp_fld + os.sep + self.id
+ self.image_list = self.parse_folder(self.folder, self.tile_id)
+
+ if len(self.image_list) > 0:
+ self.tile_id = self._tile_id(self.image_list[0])
+ self.input_dates = [self._img_date(x) for x in self.image_list]
+ self.tile_cloud_percentage = []
+
+ if input_date_interval is not None:
+ idx = [i for i in range(len(self.input_dates)) if self.input_dates[i]>=input_date_interval[0] and 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]
+
+ if not dummy_read:
+
+ if filter_by_roi is not None:
+ assert(os.path.exists(filter_by_roi))
+ idx = [i for i in range(len(self.input_dates)) if self._check_roi(self.image_list[i], filter_by_roi, roi_min_surf)]
+ self.image_list = [self.image_list[i] for i in idx]
+ self.input_dates = [self.input_dates[i] for i in idx]
+
+ if max_clouds_percentage is not None:
+ self.tile_cloud_percentage = [self._tile_cloud_percentage(x) for x in self.image_list]
+ idx = [i for i in range(len(self.input_dates)) if self.tile_cloud_percentage[i] <= max_clouds_percentage]
+ self.image_list = [self.image_list[i] for i in idx]
+ self.input_dates = [self.input_dates[i] for i in idx]
+ self.tile_cloud_percentage = [self.tile_cloud_percentage[i] for i in idx]
+
+ self.set_input_epsg()
+ self.output_epsg = self.input_epsg
+ self.output_dates = self.input_dates
+
+ for img in self.image_list:
+ for p in self.PTRN_ful:
+ ifn, ofn = self.get_file(img, p)
+ self.append(to_otb_pipeline(ifn), ofn, self.REF_TYPE, 'out', is_output=True)
+ for p in self.PTRN_msk:
+ ifn, ofn = self.get_file(img, p)
+ self.append(to_otb_pipeline(ifn), ofn, self.MSK_TYPE, 'out', is_output=True)
+
+ else:
+ warnings.warn('Empty pipeline. Need to set preprocessed inputs?')
+
+ def __del__(self):
+ if os.path.exists(self.temp_fld):
+ shutil.rmtree(self.temp_fld)
+
+ def reset(self):
+ self.pipe = []
+ self.files = []
+ self.types = []
+ self.out_p = []
+ self.out_idx = []
+
+ for img in self.image_list:
+ for p in self.PTRN_ful:
+ ifn, ofn = self.get_file(img, p)
+ self.append(to_otb_pipeline(ifn), ofn, self.REF_TYPE, 'out', is_output=True)
+ for p in self.PTRN_msk:
+ ifn, ofn = self.get_file(img, p)
+ self.append(to_otb_pipeline(ifn), ofn, self.MSK_TYPE, 'out', is_output=True)
+
+ def get_file(self, img, ptrn):
+ in_fn, out_fn = None, None
+ if os.path.isdir(img):
+ in_fn = glob.glob(os.path.join(img, ptrn))[0]
+ out_fn = in_fn.replace(self.folder, '').lstrip(os.sep)
+ elif os.path.splitext(img)[-1] == '.zip':
+ z = zipfile.ZipFile(img)
+ out_fn = [x for x in z.namelist() if ptrn.split(os.sep)[-1].replace('*','') in x][0]
+ in_fn = os.sep + 'vsizip' + os.sep + os.path.abspath(img) + os.sep + out_fn
+ return in_fn, os.path.splitext(self.tile_id + os.sep + out_fn)[0] + '.tif'
+
+ def parse_folder(self, fld, tile):
+ 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._tile_id(x) == tile]
+ zip_list = [os.path.abspath(x) for x in glob.glob(os.path.join(fld, self.PTRN_dir))
+ if os.path.splitext(x)[-1] == '.zip' and self._tile_id(x) == tile]
+ im_dict = {}
+ for i in img_list:
+ im_dict[self._img_id(i)] = i
+ for i in zip_list:
+ if self._img_id(i) not in im_dict.keys():
+ im_dict[self._img_id(i)] = i
+ return sorted(im_dict.values(), key=lambda x: self._img_id(x))
+
+ def get_coverage(self):
+ proc_idx = self.out_idx.copy()
+ proc_idx = proc_idx[1::2]
+ cc = []
+ for t in proc_idx:
+ arr = self.pipe[t].GetImageAsNumpyArray('out')
+ cc.append(np.sum(arr)/(arr.shape[0]*arr.shape[1]))
+ return cc
+
+ def set_input_epsg(self):
+ f, _ = self.get_file(self.image_list[0], self.PTRN_20m[0])
+ ds = gdal.Open(f)
+ self.input_epsg = osr.SpatialReference(wkt=ds.GetProjection()).GetAuthorityCode('PROJCS')
+ ds = None
+
+ 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 clip(self, roi):
+ assert(os.path.exists(roi))
+ proc_idx = self.out_idx.copy()
+ 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.vect', roi)
+ er.Execute()
+ fn = self.files[t]
+ ty = self.types[t]
+ self.append(er, fn, ty, 'out', is_output=True)
+
+ def preprocess(self):
+ proc_idx = self.out_idx.copy()
+ self.out_idx = []
+ num_of_mask_inputs = len(self.PTRN_msk)
+ for t in proc_idx[::10+num_of_mask_inputs]:
+ idx_to_stack = []
+ for k in range(t,t+10):
+ if k % (10 + num_of_mask_inputs) in [3, 4, 5, 7, 8, 9]:
+ si = otb.Registry.CreateApplication('Superimpose')
+ si.SetParameterInputImage('inm', self.pipe[k].GetParameterOutputImage('out'))
+ si.SetParameterInputImage('inr', self.pipe[t].GetParameterOutputImage('out'))
+ si.Execute()
+ self.append(si)
+ cr = otb.Registry.CreateApplication('BandMath')
+ cr.AddImageToParameterInputImageList('il', self.pipe[-1].GetParameterOutputImage('out'))
+ cr.SetParameterString('exp', '(im1b1 < 0 && im1b1 !=' + str(self.NDT) + ') ? 0 : im1b1')
+ cr.Execute()
+ self.append(cr)
+ idx_to_stack.append(len(self.pipe)-1)
+ else:
+ idx_to_stack.append(k)
+
+ cct = otb.Registry.CreateApplication('ConcatenateImages')
+ [cct.AddImageToParameterInputImageList('il', self.pipe[j].GetParameterOutputImage('out')) for j in idx_to_stack]
+ cct.Execute()
+ fn = self.files[t].replace(self.B2_name, 'STACK')
+ ty = self.REF_TYPE
+ self.append(cct, fn, ty, 'out', is_output=True)
+
+ msks = self.pipe[t+10:t+10+num_of_mask_inputs]
+ msk_pipe = self._process_mask(msks)
+ for app in msk_pipe[:-1]:
+ self.append(app)
+ fn = self.files[t].replace(self.B2_name, 'BINARY_MASK')
+ ty = self.MSK_TYPE
+ self.append(msk_pipe[-1], fn, ty, 'out', is_output=True)
+
+ def write_src_quicklooks(self, fld, bnds = [3,2,1], scale_factor=0.2):
+ proc_idx = self.out_idx.copy()
+ fns = [fld + os.sep + self.NAME + '_' + self._img_id(x) + '.png' for x in self.image_list]
+
+ for t,m in zip(proc_idx[::2],proc_idx[1::2]):
+ rt = otb.Registry.CreateApplication('RigidTransformResample')
+ rt.SetParameterInputImage('in', self.pipe[t].GetParameterOutputImage('out'))
+ rt.SetParameterFloat('transform.type.id.scalex', scale_factor)
+ rt.SetParameterFloat('transform.type.id.scaley', scale_factor)
+ rt.Execute()
+
+ rtm = otb.Registry.CreateApplication('RigidTransformResample')
+ rtm.SetParameterInputImage('in', self.pipe[m].GetParameterOutputImage('out'))
+ rtm.SetParameterFloat('transform.type.id.scalex', scale_factor)
+ rtm.SetParameterFloat('transform.type.id.scaley', scale_factor)
+ rtm.SetParameterString('interpolator', 'nn')
+ rtm.Execute()
+
+ bmm = otb.Registry.CreateApplication('BandMath')
+ bmm.AddImageToParameterInputImageList('il', rtm.GetParameterOutputImage('out'))
+ bmm.SetParameterString('exp', '1 - im1b1')
+ if not os.path.exists(self.temp_fld):
+ os.makedirs(self.temp_fld)
+ bmm.SetParameterString('out', self.temp_fld + os.sep + 'msk.tif')
+ bmm.ExecuteAndWriteOutput()
+
+ dc = otb.Registry.CreateApplication('DynamicConvert')
+ dc.SetParameterInputImage('in', rt.GetParameterOutputImage('out'))
+ dc.SetParameterString('mask', self.temp_fld + os.sep + 'msk.tif')
+ dc.SetParameterString('channels', 'rgb')
+ dc.SetParameterInt('channels.rgb.red', bnds[0])
+ dc.SetParameterInt('channels.rgb.green', bnds[1])
+ dc.SetParameterInt('channels.rgb.blue', bnds[2])
+ dc.SetParameterString('quantile.high', '2')
+ dc.SetParameterString('quantile.low', '2')
+ dc.SetParameterString('outmin', '1')
+ dc.SetParameterString('outmax', '255')
+
+ dc.SetParameterString('out', fns.pop(0))
+ dc.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint8)
+ dc.ExecuteAndWriteOutput()
+
+ def reproject(self, epsg):
+ if epsg != self.input_epsg:
+ self.output_epsg = epsg
+ proc_idx = self.out_idx.copy()
+ self.out_idx = []
+ i = 0
+ for t in proc_idx:
+ rp = otb.Registry.CreateApplication('OrthoRectification')
+ rp.SetParameterInputImage('io.in', self.pipe[t].GetParameterOutputImage('out'))
+ rp.SetParameterString('map', 'epsg')
+ rp.SetParameterString('map.epsg.code', self.output_epsg)
+ rp.SetParameterString('opt.gridspacing', '40')
+ if i % 2 == 0:
+ rp.SetParameterString('outputs.default', str(self.NDT))
+ fn = self.files[t].replace('_STACK.tif', '_STACK_' + self.output_epsg + '.tif')
+ ty = self.REF_TYPE
+ else:
+ rp.SetParameterString('interpolator', 'nn')
+ fn = self.files[t].replace('_BINARY_MASK.tif', '_BINARY_MASK_' + self.output_epsg + '.tif')
+ ty = self.MSK_TYPE
+ rp.Execute()
+ self.append(rp, fn, ty, 'io.out', is_output=True)
+ i += 1
+
+ def coregister(self, ref_img, ref_bnd, tgt_bnd, out_fld=None):
+ # PointMatchCoregistration is a pipeline-breaking app
+ # Need to write outputs, compute coregistered and write again
+ # Then update pipeline
+ to_del = self.write_outputs(self.temp_fld, update_pipe=True)
+ proc_idx = self.out_idx.copy()
+ self.out_idx = []
+ i = 0
+ curr_json = ''
+ for t in proc_idx:
+ cr = otb.Registry.CreateApplication('PointMatchCoregistration')
+ cr.SetParameterInputImage('in', self.pipe[t].GetParameterOutputImage('out'))
+ cr.SetParameterInt('fallback', 1)
+ if (i % 2) == 0:
+ curr_json = os.path.splitext(self.files[t])[0] + '.json'
+ curr_json = os.path.join(self.temp_fld, curr_json)
+ cr.SetParameterInt('band', tgt_bnd)
+ cr.SetParameterString('inref', ref_img)
+ cr.SetParameterInt('bandref', ref_bnd)
+ cr.SetParameterString('outjson', curr_json)
+ fn = self.files[t].replace('_STACK.tif', '_STACK_COREG.tif')
+ ty = self.REF_TYPE
+ else:
+ cr.SetParameterString('inmodel', curr_json)
+ cr.SetParameterString('interpolator', 'nn')
+ fn = self.files[t].replace('_BINARY_MASK.tif', '_BINARY_MASK_COREG.tif')
+ ty = self.MSK_TYPE
+ self.append(cr, fn, ty, 'out', is_output=True)
+ i += 1
+ if out_fld is None:
+ out_fld = self.temp_fld
+ self.write_outputs(out_fld, update_pipe=True, flag_nodata=True)
+
+ def parse_dates(self, fn):
+ with open(fn) as f:
+ return [l for l in f.read().splitlines() if l]
+
+ def gapfill(self, output_dates=None, on_disk=False):
+ #assert(os.path.exists(output_dates))
+ proc_idx = self.out_idx.copy()
+ self.out_idx = []
+
+ stk = otb.Registry.CreateApplication('ConcatenateImages')
+ [stk.AddImageToParameterInputImageList('il', self.pipe[x].GetParameterOutputImage(self.out_p[x])) for x in proc_idx[::2]]
+ stk.Execute()
+ self.append(stk)
+
+ mst = otb.Registry.CreateApplication('ConcatenateImages')
+ [mst.AddImageToParameterInputImageList('il', self.pipe[x].GetParameterOutputImage(self.out_p[x])) for x in proc_idx[1::2]]
+ mst.Execute()
+ self.append(mst)
+
+ with open(os.path.join(self.folder, self.tile_id + '_indates.txt'), 'w') as df:
+ [df.write(x + '\n') for x in self.input_dates]
+
+ if output_dates is not None:
+ od = self.parse_dates(output_dates)
+ self.output_dates = od
+ else:
+ output_dates = os.path.join(self.folder, self.tile_id + '_indates.txt')
+ od = self.input_dates
+
+ self.output_dates = od
+
+ gf = otb.Registry.CreateApplication('ImageTimeSeriesGapFilling')
+ gf.SetParameterInputImage('in', self.pipe[-2].GetParameterOutputImage('out'))
+ gf.SetParameterInputImage('mask', self.pipe[-1].GetParameterOutputImage('out'))
+ gf.SetParameterInt('comp', 10)
+ gf.SetParameterString('it', 'linear')
+ gf.SetParameterString('id', os.path.join(self.folder, self.tile_id + '_indates.txt'))
+ gf.SetParameterString('od', output_dates)
+ if not on_disk:
+ gf.Execute()
+ gf_out = gf
+ else:
+ gf_fn = self.temp_fld + os.sep + 'SENTINEL2_' + self.tile_id + '_GAPFILLED_FULL.tif'
+ if not os.path.exists(gf_fn):
+ if not os.path.exists(self.temp_fld):
+ os.makedirs(self.temp_fld)
+ gf.SetParameterString("out", gf_fn)
+ gf.ExecuteAndWriteOutput()
+ gf_out = to_otb_pipeline(gf_fn)
+ self.append(gf_out)
+
+ t = 1
+ for d in od:
+ ch_list = ['Channel%d' % i for i in range(t,t+10)]
+ t += 10
+ er = otb.Registry.CreateApplication('ExtractROI')
+ er.SetParameterInputImage('in', gf_out.GetParameterOutputImage('out'))
+ er.UpdateParameters()
+ er.SetParameterStringList('cl', ch_list)
+ er.Execute()
+ dn = 'SENTINEL2_' + self.tile_id + '_GAPFILL_' + d
+ fn = self.tile_id + os.sep + dn + os.sep + dn + '_STACK.tif'
+ ty = self.REF_TYPE
+ self.append(er, fn, ty, 'out', is_output=True)
+
+ def generate_feature_stack(self, feat_list=None):
+ proc_idx = self.out_idx.copy()
+ self.out_idx = []
+
+ exp_list = self.FEAT_exp.values()
+ stack_name = 'FEAT'
+ if feat_list is not None:
+ exp_list = [self.FEAT_exp[x] for x in feat_list]
+ stack_name = '_'.join(feat_list)
+
+ for t in proc_idx:
+ bm = otb.Registry.CreateApplication('BandMathX')
+ bm.AddImageToParameterInputImageList('il', self.pipe[t].GetParameterOutputImage('out'))
+ bm.SetParameterString('exp', '{' + ';'.join(exp_list) + '}')
+ bm.Execute()
+ fn = self.files[t].replace('_STACK.tif', '_' + stack_name + '.tif')
+ ty = otb.ImagePixelType_float
+ self.append(bm, fn, ty, 'out', is_output=True)
+
+ return stack_name
+
+ def generate_time_series(self, feat_list):
+ proc_idx = self.out_idx.copy()
+ self.out_idx = []
+
+ for feat in feat_list:
+ bm = otb.Registry.CreateApplication('BandMathX')
+ i = 1
+ expr = []
+ for t in proc_idx:
+ expr.append(self.FEAT_exp[feat].replace('im1', 'im%d' % i))
+ bm.AddImageToParameterInputImageList('il', self.pipe[t].GetParameterOutputImage('out'))
+ i += 1
+ bm.SetParameterString('exp', '{' + ';'.join(expr) + '}')
+ bm.Execute()
+ fn = self.tile_id + os.sep + 'SENTINEL2_' + self.tile_id + '_GAPFILL_' + feat + '.tif'
+ ty = otb.ImagePixelType_float
+ self.append(bm, fn, ty, 'out', is_output=True)
+
+ def set_preprocess_output(self, fld):
+ lst = self.parse_folder(os.path.join(fld, self.tile_id), self.tile_id)
+ if len(lst) == len(self.image_list):
+ self.out_idx = []
+ for img in lst:
+ f = glob.glob(os.path.join(img, '*STACK.tif'))[0]
+ self.append(to_otb_pipeline(f), f, self.REF_TYPE, 'out', is_output=True)
+ f = glob.glob(os.path.join(img, '*BINARY_MASK.tif'))[0]
+ self.append(to_otb_pipeline(f), f, self.MSK_TYPE, 'out', is_output=True)
+ else:
+ warnings.warn("No matching files found as preprocess output. No modification to pipeline.")
+
+ def set_gapfilled_output(self, file_pattern):
+ lst = sorted(glob.glob(file_pattern))
+ self.out_idx = []
+ if len(lst) > 0:
+ for f in lst:
+ pf = os.path.join(self.tile_id, os.path.basename(f))
+ self.append(to_otb_pipeline(f), pf, self.REF_TYPE, 'out', is_output=True)
+ else:
+ warnings.warn("No matching files found as preprocess output. No modification to pipeline.")
+
+ def write_outputs(self, fld, roi=None, update_pipe=False, compress=False, flag_nodata=False):
+ out = []
+ proc_idx = self.out_idx.copy()
+ if roi is not None:
+ out_idx_bck = self.out_idx.copy()
+ pipe_length = len(self.pipe)
+ self.clip(roi)
+ if update_pipe:
+ assert roi is None, 'Cannot set output files as pipe input over a ROI, use clip function instead.'
+ 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)
+ if roi is not None:
+ self.out_idx = out_idx_bck
+ self.pipe = self.pipe[:pipe_length]
+ self.files = self.files[:pipe_length]
+ self.types = self.types[:pipe_length]
+ self.out_p = self.out_p[:pipe_length]
+
+ if flag_nodata:
+ if isinstance(flag_nodata, bool):
+ val = self.NDT
+ elif isinstance(flag_nodata, int) or isinstance(flag_nodata,float):
+ val = flag_nodata
+ for f in out:
+ ds = gdal.Open(f, 1)
+ for i in range(ds.RasterCount):
+ ds.GetRasterBand(i+1).SetNoDataValue(val)
+ ds = None
+
+ return out
+
+
+class S2TheiaPipeline:
+
+ S2TilePipeline = S2TheiaTilePipeline
+ _check = S2TilePipeline._check
+ _tile_id = S2TilePipeline._tile_id
+
+ def __init__(self, fld, temp_fld='/tmp', input_date_interval=None, max_clouds_percentage=None):
+ self.folder = fld
+ self.temp_fld = temp_fld
+ self.input_date_interval = input_date_interval
+ self.max_clouds_percentage = max_clouds_percentage
+ self.tile_list = set()
+ img_list = [os.path.abspath(x) for x in glob.glob(os.path.join(self.folder, self.S2TilePipeline.PTRN_dir))
+ if os.path.isdir(x) or os.path.splitext(x)[-1] == '.zip']
+ [self.tile_list.add(self._tile_id(x)) for x in img_list if self._check(x)]
+ self.tiles = [self.S2TilePipeline(fld, t, self.temp_fld, self.input_date_interval, self.max_clouds_percentage) for t in self.tile_list]
+ self.output_dates = None
+ self.roi = None
+ self.output_epsg = self.tiles[0].input_epsg
+
+ def __del__(self):
+ for x in self.tiles:
+ del x
+
+ def set_output_dates_by_file(self, od):
+ self.output_dates = od
+
+ def set_output_dates(self, start, end, step=10):
+ start_date = datetime.datetime(int(start[0:4]), int(start[4:6]), int(start[6:8]))
+ end_date = datetime.datetime(int(end[0:4]), int(end[4:6]), int(end[6:8]))
+ d, st = start_date, datetime.timedelta(step)
+ tmstp = str(datetime.datetime.timestamp(datetime.datetime.now())).replace('.', '')
+ ofn = self.folder + os.sep + 's2ppl_' + tmstp + '_output_dates.txt'
+ with open(ofn, 'w') as f:
+ while d < end_date:
+ f.write(d.strftime('%Y%m%d')+'\n')
+ d += st
+ self.output_dates = ofn
+
+ def set_roi(self, roi):
+ self.roi = roi
+
+ def set_output_epsg(self, epsg):
+ self.output_epsg = epsg
+
+ def extract_feature_set(self, out_fld, feat_list=None, mosaicking=None, store_gapfill=False):
+ out = []
+ stack_name = ''
+ if self.output_dates is not None:
+ for t in self.tiles:
+ if self.roi is not None:
+ t.clip(self.roi)
+ t.preprocess()
+ t.reproject(self.output_epsg)
+ t.gapfill(self.output_dates, store_gapfill)
+ stack_name = t.generate_feature_stack(feat_list)
+ out.append(t.write_outputs(out_fld))
+ t.reset()
+ if mosaicking == 'vrt':
+ out_mos = []
+ vrtopt = gdal.BuildVRTOptions()
+ for i in range(len(self.tiles[0].output_dates)):
+ fn = out_fld + os.sep + 'SENTINEL2_MOSAIC_GAPFILL_' + self.tiles[0].output_dates[i] + '_' + stack_name + '.vrt'
+ to_mosaic = [x[i] for x in out]
+ gdal.BuildVRT(fn, to_mosaic, options=vrtopt)
+ out_mos.append(fn)
+ return out_mos
+ return out
+
+ def extract_time_series(self, out_fld, feat_list, mosaicking=None, store_gapfill=False):
+ out = []
+ if self.output_dates is not None:
+ for t in self.tiles:
+ if self.roi is not None:
+ t.clip(self.roi)
+ t.preprocess()
+ t.reproject(self.output_epsg)
+ t.gapfill(self.output_dates, store_gapfill)
+ t.generate_time_series(feat_list)
+ out.append(t.write_outputs(out_fld))
+ t.reset()
+ if mosaicking == 'vrt':
+ out_mos = []
+ vrtopt = gdal.BuildVRTOptions()
+ for i in range(len(feat_list)):
+ fn = out_fld + os.sep + 'SENTINEL2_MOSAIC_GAPFILL_' + feat_list[i] + '.vrt'
+ to_mosaic = [x[i] for x in out]
+ gdal.BuildVRT(fn, to_mosaic, options=vrtopt)
+ out_mos.append(fn)
+ return out_mos
+ return out
+
+
+"""
+def single_tile_clip(ppl : S2Pipeline, roi, write_output=False, s2out=None):
+
+ if write_output and s2out is None:
+ sys.exit("Please provide output folder (s2out=) when write_output is true")
+
+ if s2out is not None:
+ s2out = os.path.abspath(s2out)
+
+ out_list = []
+ for l in ppl.pipe:
+ cur_list = []
+ for f in l:
+ er = otb.Registry.CreateApplication('ExtractROI')
+ er.SetParameterString('in', f)
+ er.SetParameterString('mode', 'fit')
+ er.SetParameterString('mode.fit.vect', roi)
+ if write_output:
+ dst_file = os.path.join(s2out, f.replace(ppl.s2_folder, ''))
+ if not os.path.exists(os.path.dirname(dst_file)):
+ os.makedirs(os.path.dirname(dst_file))
+ typ = MSK_TYPE
+ if PTRN_ref in os.path.basename(dst_file):
+ typ = REF_TYPE
+ er.SetParameterString('out', dst_file)
+ er.SetParameterOutputImagePixelType('out', typ)
+ er.ExecuteAndWriteOutput()
+ cur_list.append(to_otb_pipeline(dst_file))
+ else:
+ er.Execute()
+ cur_list.append(er)
+ out_list.append(cur_list)
+
+ ppl.set_pipeline(out_list)
+
+ return
+
+def single_tile_preprocess(ppl : S2Pipeline, write_output=False, output_file=None):
+
+ ppl, dl = single_tile_scan(s2fld)
+
+ begin_index = len(ppl)
+
+ for d in ppl[:begin_index]:
+ for i in range(10):
+ if i in [3, 4, 5, 7, 8, 9]:
+ ppl.append(otb.Registry.CreateApplication('Superimpose'))
+ ppl[-1].SetParameterInputImage('inm', d[i].GetParameterOutputImage('out'))
+ ppl[-1].SetParameterInputImage('inr', d[0].GetParameterOutputImage('out'))
+ ppl[-1].Execute()
+ else:
+ ppl.append(d[i])
+
+ ppl.append(otb.Registry.CreateApplication('ConcatenateImages'))
+ [ppl[-1].AddImageToParameterInputImageList('il', x.GetParameterOutputImage('out')) for x in ppl[begin_index:-1]]
+ ppl[-1].Execute()
+
+ stack_index = len(ppl) - 1
+
+ for d in ppl[:begin_index]:
+ ppl.append(otb.Registry.CreateApplication('BandMath'))
+ [ppl[-1].AddImageToParameterInputImageList('il', x.GetParameterOutputImage('out')) for x in d[10:]]
+ ppl[-1].SetParameterString('exp', 'im1b1!=0 || im2b1!=0 || im3b1!=0')
+ ppl[-1].Execute()
+
+ ppl.append(otb.Registry.CreateApplication('ConcatenateImages'))
+ [ppl[-1].AddImageToParameterInputImageList('il', x.GetParameterOutputImage('out')) for x in ppl[stack_index+1:-1]]
+ ppl[-1].Execute()
+
+ mask_index = len(ppl) - 1
+
+ with open(os.path.join(s2fld, 'indates.txt'), 'w') as df:
+ [df.write(x + '\n') for x in dl]
+
+ ppl.append(otb.Registry.CreateApplication('ImageTimeSeriesGapFilling'))
+ ppl[-1].SetParameterInputImage('in', ppl[stack_index].GetParameterOutputImage('out'))
+ ppl[-1].SetParameterInputImage('mask', ppl[mask_index].GetParameterOutputImage('out'))
+ ppl[-1].SetParameterInt('comp', 10)
+ ppl[-1].SetParameterString('it', 'linear')
+ ppl[-1].SetParameterString('id', os.path.join(s2fld, 'indates.txt'))
+ ppl[-1].SetParameterString('od', output_dates)
+
+ if write_output is False:
+ ppl[-1].Execute()
+ return ppl
+ else:
+ if output_file is None:
+ sys.exit("Please provide an output file name when write_output is True.")
+ ppl[-1].SetParameterString('out', output_file)
+ ppl[-1].SetParameterOutputImagePixelType('out', REF_TYPE)
+ ppl[-1].ExecuteAndWriteOutput()
+ return output_file
+
+def single_tile_gapfilling()
+
+def single_tile_features(gf_stack, output_dates, write_output=False, output_folder=None):
+
+ ppl = to_otb_pipeline(gf_stack)
+ n_dates = int(ppl[-1].GetImageNbBands('out') / 10)
+ print(n_dates)
+
+ with open(output_dates,'r') as df:
+ dates = df.read().splitlines()
+
+ if n_dates != len(dates):
+ sys.exit('Number of dates in stack does not correspond to the one in the text file.')
+
+ gf_index = len(ppl) - 1
+
+ # Extract dates from gapfilling stacks
+ for d in range(n_dates):
+ ppl.append(otb.Registry.CreateApplication('ExtractROI'))
+ ppl[-1].SetParameterInputImage('in', ppl[gf_index].GetParameterOutputImage('out'))
+ ppl[-1].UpdateParameters()
+ ppl[-1].SetParameterStringList('cl', ['Channel%d' % i for i in range(10*d+1, 10*d+11)])
+ ppl[-1].Execute()
+
+ out_files = []
+
+ # Compute per-date features
+ for d in range(n_dates):
+ ppl.append(otb.Registry.CreateApplication('BandMathX'))
+ ppl[-1].AddImageToParameterInputImageList('il', ppl[gf_index+1+d].GetParameterOutputImage('out'))
+ ppl[-1].SetParameterString('exp', '{' + ';'.join(FEAT_exp) + '}')
+ if write_output is True:
+ print(True)
+ else:
+ ppl[-1].Execute()
+
+ if write_output is True:
+ return out_files, n_dates
+ else:
+ return ppl, n_dates
+
+
+if __name__ == "__main__":
+ single_tile_clip(sys.argv[1],sys.argv[2],sys.argv[3])
+"""
\ No newline at end of file
diff --git a/VHR/vhrbase.py b/VHR/vhrbase.py
new file mode 100644
index 0000000000000000000000000000000000000000..175bf091840d81a78a10ad69e3223661a97b4092
--- /dev/null
+++ b/VHR/vhrbase.py
@@ -0,0 +1,189 @@
+import otbApplication as otb
+from otb_numpy_proc import to_otb_pipeline
+import os
+import glob
+from pathlib import Path
+import re
+import ogr
+
+class SPOT67RasterPipeline:
+ # BEGIN SPOT6/7 VHR PROTOTYPE
+ MS_FOLDER_PATTERN = 'IMG_SPOT*_MS_*'
+ PAN_FOLDER_PATTERN = 'IMG_SPOT*_P_*'
+ TILE_PATTERN = '*_R*C*.TIF'
+ NDT = 0.0
+ REF_TYPE = otb.ImagePixelType_uint16
+ TOA_TYPE = otb.ImagePixelType_float
+ GRID_SPC = [6, 24]
+ RS_LABEL = 'SEN'
+ MS_LABEL = 'MS'
+ OPT_IN = '?&skipcarto=true'
+
+ @classmethod
+ def _find_tiles(cls, x):
+ lst = sorted(glob.glob(os.path.join(x, cls.TILE_PATTERN)))
+ res = re.search('R[0-9]C[0-9]', lst[-1]).group()
+ return lst, int(res[1]), int(res[3]), res
+
+ def __init__(self, pan_root, ms_root=None):
+ self.out_idx = []
+ self.pipe = []
+ self.files = []
+ self.types = []
+ self.out_p = []
+
+ self.ms_folder = None
+ self.pan_folder = None
+ self.pipe = []
+
+ if ms_root is None:
+ ms_root = pan_root
+
+ for p in Path(pan_root).rglob(self.PAN_FOLDER_PATTERN):
+ if os.path.isdir(p):
+ self.pan_folder = str(p)
+ break
+
+ for p in Path(ms_root).rglob(self.MS_FOLDER_PATTERN):
+ if os.path.isdir(p):
+ self.ms_folder = str(p)
+ break
+
+ if self.ms_folder is None or self.pan_folder is None:
+ raise ValueError('No MS and/or PAN folder found')
+
+ self.pan_list, r_pan, c_pan, res = self._find_tiles(self.pan_folder)
+ tile_p = otb.Registry.CreateApplication('TileFusion')
+ lst = [x + self.OPT_IN for x in self.pan_list]
+ tile_p.SetParameterStringList('il', lst)
+ tile_p.SetParameterInt('rows', r_pan)
+ tile_p.SetParameterInt('cols', c_pan)
+ tile_p.Execute()
+ fn = os.path.basename(self.pan_list[-1]).replace(res, 'FULL')
+ ty = self.REF_TYPE
+ self.append(tile_p, fn, ty, 'out', is_output=True)
+
+ self.ms_list, r_ms, c_ms, res = self._find_tiles(self.ms_folder)
+ tile_ms = otb.Registry.CreateApplication('TileFusion')
+ lst = [x + self.OPT_IN for x in self.ms_list]
+ tile_ms.SetParameterStringList('il', lst)
+ tile_ms.SetParameterInt('rows', r_ms)
+ tile_ms.SetParameterInt('cols', c_ms)
+ tile_ms.Execute()
+ fn = os.path.basename(self.ms_list[-1]).replace(res, 'FULL')
+ ty = self.REF_TYPE
+ self.append(tile_ms, fn, ty, 'out', is_output=True)
+
+ def to_toa(self, clamp=True):
+ proc_idx = self.out_idx.copy()
+ self.out_idx = []
+
+ for p in proc_idx:
+ fn = self.files[p].replace('_FULL', '_FULL_TOA')
+ toa = otb.Registry.CreateApplication('OpticalCalibration')
+ toa.SetParameterInputImage('in', self.pipe[p].GetParameterOutputImage(self.out_p[p]))
+ toa.SetParameterInt('clamp', int(clamp))
+ toa.Execute()
+ self.append(toa)
+ bm = otb.Registry.CreateApplication('BandMathX')
+ bm.AddImageToParameterInputImageList('il', self.pipe[-1].GetParameterOutputImage('out'))
+ bm.SetParameterString('exp', '10000*im1')
+ bm.Execute()
+ ty = self.REF_TYPE
+ self.append(bm, fn, ty, 'out', is_output=True)
+ return
+
+ def orthorectify(self, dem_fld, geoid = None):
+ assert (os.path.isdir(dem_fld))
+ proc_idx = self.out_idx.copy()
+ self.out_idx = []
+ gspc_idx = 0
+ 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', self.GRID_SPC[gspc_idx])
+ ortho.Execute()
+ fn = self.files[t].replace(self.RS_LABEL, 'ORTHO')
+ ty = self.REF_TYPE
+ self.append(ortho, fn, ty, 'io.out', is_output=True)
+ gspc_idx += 1
+
+ def pansharp(self):
+ proc_idx = self.out_idx.copy()
+ self.out_idx = []
+ assert(len(proc_idx) == 2)
+ btps = otb.Registry.CreateApplication("BundleToPerfectSensor")
+ btps.SetParameterInputImage('inp', self.pipe[proc_idx[0]].GetParameterOutputImage(self.out_p[proc_idx[0]]))
+ btps.SetParameterInputImage('inxs', self.pipe[proc_idx[1]].GetParameterOutputImage(self.out_p[proc_idx[1]]))
+ btps.SetParameterString('method', 'bayes')
+ btps.Execute()
+ fn = self.files[proc_idx[1]].replace(self.MS_LABEL, 'PS')
+ ty = self.REF_TYPE
+ self.append(btps, fn, ty, 'out', is_output=True)
+
+ def clip(self, roi, buffer=0):
+ assert(os.path.exists(roi))
+ proc_idx = self.out_idx.copy()
+ 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', 'extent')
+ ds = ogr.Open(roi)
+ ly = ds.GetLayer(0)
+ extent = ly.GetExtent()
+ extent = (extent[0] - buffer, extent[1] + buffer, extent[2] - buffer, extent[3] + buffer)
+ if ly.GetSpatialRef().GetAuthorityCode('PROJCS') == '4326':
+ er.SetParameterString('mode.extent.unit', 'lonlat')
+ else:
+ er.SetParameterString('mode.extent.unit', 'phy')
+ er.SetParameterFloat('mode.extent.ulx', extent[0])
+ er.SetParameterFloat('mode.extent.uly', extent[2])
+ er.SetParameterFloat('mode.extent.lrx', extent[1])
+ er.SetParameterFloat('mode.extent.lry', extent[3])
+ ly, ds = None, None
+ er.Execute()
+ fn = self.files[t]
+ ty = self.types[t]
+ self.append(er, fn, ty, 'out', is_output=True)
+
+ 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 write_outputs(self, fld, roi=None, update_pipe=False, compress=False):
+ out = []
+ out_idx_bck = self.out_idx.copy()
+ if roi is not None:
+ pipe_length = len(self.pipe)
+ self.clip(roi)
+ if update_pipe:
+ assert roi is None, 'Cannot set output files as pipe input over a ROI.'
+ self.out_idx = []
+ for t in out_idx_bck:
+ 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)
+ if roi is not None:
+ self.out_idx = out_idx_bck
+ self.pipe = self.pipe[:pipe_length]
+ self.files = self.files[:pipe_length]
+ self.types = self.types[:pipe_length]
+ self.out_p = self.out_p[:pipe_length]
+ return out
\ No newline at end of file