diff --git a/Common/otb_numpy_proc.py b/Common/otb_numpy_proc.py index 00b6bdac2d73ab3b389ed743849071ba189f877d..81beba3fc6e2a99c083ca734e037e08f527bcafd 100644 --- a/Common/otb_numpy_proc.py +++ b/Common/otb_numpy_proc.py @@ -1,4 +1,3 @@ -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)) diff --git a/OBIA/OBIABase.py b/OBIA/OBIABase.py index 9f9207992e37a70d122711dd48cb131a8d434eb4..b2ab8afa9bf3817014217dcaef6f6e3cbb975dbc 100644 --- a/OBIA/OBIABase.py +++ b/OBIA/OBIABase.py @@ -1,6 +1,5 @@ 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' diff --git a/OBIA/segmentation.py b/OBIA/segmentation.py index 5be2e7e19f89e60f0e25d633dc0a829b83d14e92..752455b590d0d49f33c021f2390efcd5b9cda611 100644 --- a/OBIA/segmentation.py +++ b/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]), diff --git a/TimeSeries/s2theia.py b/TimeSeries/s2theia.py index 9c6857aad0eed5cfab20d8596e7c10a16c62c7b8..78e847491f12a8f9427e9fc568b25296567e32af 100644 --- a/TimeSeries/s2theia.py +++ b/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)): diff --git a/moringa.py b/moringa.py index ca764aade1256040f2eb25efa0e833b4456c9598..e20905f195213df9eb009561e9ce999b0b3c653c 100755 --- a/moringa.py +++ b/moringa.py @@ -1,7 +1,9 @@ +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,