ENH: added preprocess_s2 workflow into main script (plus minor changes).

Common/otb_numpy_proc.py

-import otbApplication as otb
 import numpy as np
 import sys
 from time import sleep
@@ -6,6 +5,7 @@ import os, psutil
 from functools import partial
 import string, random
 from osgeo.gdal import BuildVRT
+import otbApplication as otb
 
 def randomword(length):
     return ''.join(random.choice(string.ascii_lowercase) for i in range(length))

OBIA/OBIABase.py

 import os.path
 import sys
-import otbApplication as otb
 import rasterio as rio
 from rasterio.features import shapes
 import numpy as np
@@ -9,6 +8,7 @@ import pandas as pd
 from enum import Enum
 from Common.otb_numpy_proc import to_otb_pipeline
 from skimage.measure import regionprops, label
+import otbApplication as otb
 
 class OStats(Enum):
     MEAN = 'mean'

OBIA/segmentation.py

@@ -45,7 +45,7 @@ def lsgrm_process_tile(input_image, params : LSGRMParams, tile_width, tile_heigh
     seg_app.SetParameterInt('xtileidx', tile_idx[0])
     seg_app.SetParameterInt('ytileidx', tile_idx[1])
     seg_app.SetParameterString('out', out_graph)
-    seg_app.Execute()
+    seg_app.ExecuteAndWriteOutput()
 
     return [out_graph + '_node_{}_{}.bin'.format(tile_idx[1], tile_idx[0]),
             out_graph + '_nodeMargin_{}_{}.bin'.format(tile_idx[1], tile_idx[0]),

TimeSeries/s2theia.py

@@ -434,7 +434,10 @@ class S2TheiaTilePipeline:
         mst.Execute()
         self.append(mst)
 
-        with open(os.path.join(self.folder, self.tile_id + '_indates.txt'), 'w') as df:
+        if not os.path.exists(self.temp_fld):
+            os.makedirs(self.temp_fld)
+
+        with open(os.path.join(self.temp_fld, self.tile_id + '_indates.txt'), 'w') as df:
             [df.write(x + '\n') for x in self.input_dates]
 
         if output_dates is not None:
@@ -459,8 +462,6 @@ class S2TheiaTilePipeline:
         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)
@@ -616,7 +617,7 @@ class S2TheiaPipeline:
         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'
+        ofn = self.temp_fld + 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')
@@ -629,7 +630,8 @@ class S2TheiaPipeline:
     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):
+    def extract_feature_set(self, out_fld, feat_list=None, mosaicking=None, store_gapfill=False,
+                            coregister_to=None, coregister_to_band=1, coregister_using_band=3):
         out = []
         stack_name = ''
         if self.output_dates is not None:
@@ -638,11 +640,13 @@ class S2TheiaPipeline:
                     t.clip(self.roi)
                 t.preprocess()
                 t.reproject(self.output_epsg)
+                if coregister_to is not None:
+                    t.coregister(coregister_to, coregister_to_band, coregister_using_band, t.temp_fld)
                 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':
+            if len(self.tiles) > 1 and mosaicking == 'vrt':
                 out_mos = []
                 vrtopt = gdal.BuildVRTOptions()
                 for i in range(len(self.tiles[0].output_dates)):

moringa.py

+import os.path
 import sys
 import argparse
 import OBIA.segmentation
 import VHR.vhrbase
+from TimeSeries import s2theia
 
 def run_segmentation(img, threshold, cw, sw , out_seg,
                      n_first_iter, margin, n_proc, memory,
@@ -22,13 +24,40 @@ def preprocess_spot67(in_fld, out_fld, dem_fld, geoid_file, skip_ps, compress):
         sp.write_outputs(out_fld, compress=compress)
     return
 
+def preprocess_s2(in_fld, out_fld, output_dates_file=None, roi=None, coregister_to=None, coregister_to_band=1,
+                  coregister_using_band=3, provider='theia'):
+    S2Processor = None
+    if provider=='theia':
+        S2Processor = s2theia.S2TheiaPipeline
+    else:
+        raise ValueError("Unsupported/non-valid provider")
+
+    s2 = S2Processor(in_fld, temp_fld=out_fld)
+    if roi is not None:
+        s2.set_roi(roi)
+    if output_dates_file is not None:
+        s2.set_output_dates_by_file(output_dates_file)
+    else:
+        raise ValueError("Please provide path to a text file containing output dates.")
+
+    s2.extract_feature_set(out_fld, coregister_to=coregister_to, coregister_to_band=coregister_to_band,
+                           coregister_using_band=coregister_using_band, store_gapfill=True, mosaicking='vrt')
+    return
+
 def main(args):
     parser = argparse.ArgumentParser(prog="moringa", add_help=False)
     subpar = parser.add_subparsers(dest="cmd")
 
-    prepr = subpar.add_parser("preprocess", help="Performs basic time series preprocessing (unimplemented).",
+    prepr = subpar.add_parser("preprocess_s2", help="Performs Moringa preset time series preprocessing.",
                               formatter_class=argparse.ArgumentDefaultsHelpFormatter)
-    # Put arguments and options here
+    prepr.add_argument("in_folder", type=str, help="Path to the folder containing (de-zipped) S2 images.")
+    prepr.add_argument("out_folder", type=str, help="Path to the folder in which preprocessed stacks will be stored.")
+    prepr.add_argument("--output_dates_file", type=str, default=None, help="Path to the text file containing output dates for  temporal interpolation.")
+    prepr.add_argument("--roi", type=str, default=None, help="Path to the ROI vector file.")
+    prepr.add_argument("--coregister_to", type=str, default=None, help="Path to a reference image to which the stacks must be coregistered.")
+    prepr.add_argument("--coregister_to_band", type=int, nargs='?', default=1, help="Band of reference image used for co-registration.")
+    prepr.add_argument("--coregister_using_band", type=int, nargs='?', default=3, help="Band of reference image used for co-registration.")
+    prepr.add_argument("--provider", type=str, default='theia', help="S2 image provider. Supported: 'theia', 'theial3a', 'sen2cor', 'planetary'")
 
     segmt = subpar.add_parser("segment", help="Performs (large scale Baatz-Shape) segmentation of an input image.",
                               formatter_class=argparse.ArgumentDefaultsHelpFormatter)
@@ -61,8 +90,9 @@ def main(args):
 
     arg = parser.parse_args()
 
-    if arg.cmd == "preprocess":
-        print('Not yet implemented.')
+    if arg.cmd == "preprocess_s2":
+        preprocess_s2(arg.in_folder, arg.out_folder, output_dates_file=arg.output_dates_file, roi=arg.roi,
+                      coregister_to=arg.coregister_to, provider=arg.provider)
 
     if arg.cmd == "segment":
         run_segmentation(arg.img, arg.threshold, arg.cw, arg.sw, arg.outimg, arg.n_first_iter, arg.tile_margin,