TimeSeries/s2planetary.py

@@ -7,6 +7,10 @@ from Common.geotools import get_query_bbox
 from pyproj import Transformer as T
 from shapely.geometry import Polygon
 import tqdm
+from datetime import *
+import geopandas as gpd
+from scipy.interpolate import interp1d
+
 
 def fetch(shp, dt, output_fld, auth=None, band_list=None):
 
@@ -20,6 +24,7 @@ def fetch(shp, dt, output_fld, auth=None, band_list=None):
     if band_list is not None:
         lst = band_list
 
+    ofiles = []
     prg = tqdm.tqdm(total=len(res.item_collection()) * len(lst), desc="Fetching from Planetary")
     for item in res.items():
         with rasterio.open(item.assets['B02'].href) as ds:
@@ -34,6 +39,7 @@ def fetch(shp, dt, output_fld, auth=None, band_list=None):
 
         for a in lst:
             ofn = os.path.join(output_fld, '/'.join(item.assets[a].get_absolute_href().split('?')[0].split('/')[10:]))
+            ofiles.append(ofn)
             os.makedirs(os.path.dirname(ofn),exist_ok=True)
             with rasterio.open(item.assets[a].href) as img:
                 out_img, out_geot = rasterio.mask.mask(img,
@@ -54,7 +60,116 @@ def fetch(shp, dt, output_fld, auth=None, band_list=None):
         return S2PlaneteryPipeline(output_fld)
     else:
         warnings.warn("Queried for a non-default band list. Skipping pipeline setup.")
-        return
+        return ofiles
+
+def split_dt_monthly(dt):
+    sdt, edt = [datetime.strptime(x,'%Y-%m-%d') for x in dt.split('/')]
+    dts = [sdt]
+    cdt = sdt
+    while cdt < edt:
+        cdt += timedelta(weeks=4)
+        dts.append(cdt)
+    if dts[-1] > edt:
+        dts = dts[:-1]
+    if dts[-1] < edt:
+        dts.append(edt+timedelta(1))
+    out = []
+    for i in range(len(dts)-1):
+        out.append(datetime.strftime(dts[i], '%Y-%m-%d') + '/' + datetime.strftime(dts[i+1]-timedelta(1), '%Y-%m-%d'))
+    return out
+    
+def sample(shp, dt, band_list, field_name, temp_dir = '/tmp', strategy='all', auth = None):
+    fd = dt.split('/')[0]
+    fd = datetime.strptime(fd,'%Y-%m-%d')
+    ed = fd + timedelta(5)
+    ndt = datetime.strftime(fd, '%Y-%m-%d') + '/' + datetime.strftime(ed, '%Y-%m-%d')
+    #fn = fetch(shp, ndt, temp_dir, band_list=[band_list[0]])
+    fn = ['/work/temp/S2B_MSIL2A_20210104T103329_N0212_R108_T30PVT_20210121T134941.SAFE/GRANULE/L2A_T30PVT_A020011_20210104T104033/IMG_DATA/R10m/T30PVT_20210104T103329_B04_10m.tif',
+          '/work/temp/S2B_MSIL2A_20210104T103329_N0212_R108_T30PUT_20210105T133826.SAFE/GRANULE/L2A_T30PUT_A020011_20210104T104033/IMG_DATA/R10m/T30PUT_20210104T103329_B04_10m.tif']
+    
+    mos = otb.Registry.CreateApplication('Mosaic')
+    mos.SetParameterStringList('il', fn)
+    mos.Execute()
+    sts = otb.Registry.CreateApplication('PolygonClassStatistics')
+    sts.SetParameterInputImage('in', mos.GetParameterOutputImage('out'))
+    sts.SetParameterString('vec', shp)
+    sts.SetParameterString('field', field_name)
+    sts.SetParameterString('out', os.path.join(temp_dir, 'stats.xml'))
+    sts.ExecuteAndWriteOutput()
+
+    stats_fn = os.path.join(temp_dir, 'stats.xml')
+    samples_fn = os.path.join(temp_dir, 'samples.sqlite')
+
+    mos = otb.Registry.CreateApplication('Mosaic')
+    mos.SetParameterStringList('il', fn)
+    mos.Execute()
+    smp = otb.Registry.CreateApplication('SampleSelection')
+    smp.SetParameterInputImage('in', mos.GetParameterOutputImage('out'))
+    smp.SetParameterString('vec', shp)
+    smp.SetParameterString('instats', stats_fn)
+    smp.SetParameterString('field', field_name)
+    if isinstance(strategy, str):
+        smp.SetParameterString('strategy', strategy)
+    elif isinstance(strategy, int):
+        smp.SetParameterString('strategy','constant')
+        smp.SetParameterInt('strategy.constant.nb', strategy)
+    smp.SetParameterString('out', samples_fn)
+    smp.ExecuteAndWriteOutput()
+
+    ds = gpd.read_file(samples_fn)
+    coords = list(zip(list(ds.geometry.x),list(ds.geometry.y)))
+    bbox, _, _ = get_query_bbox(shp, return_all=True)
+    if auth is not None:
+        os.environ['PC_SDK_SUBSCRIPTION_KEY'] = auth
+    
+    all_arr = []
+    dates = []
+    for dti in split_dt_monthly(dt):
+
+        api = Client.open('https://planetarycomputer.microsoft.com/api/stac/v1', modifier=PC.sign_inplace)
+        res = api.search(collections="sentinel-2-l2a", bbox=bbox, datetime=dti, sortby='datetime')
+        found = res.item_collection()
+
+        if 'SCL' not in band_list:
+            band_list.append('SCL')
+        arr = {x : np.empty((ds.shape[0],len(found))) for x in band_list}
+        i = 0
+        prg = tqdm.tqdm(total=len(found), desc="Sampling from Planetary on period {}".format(dti))
+        for item in found:
+            dates.append(item.datetime)
+            for band in band_list:
+                with rasterio.open(item.assets[band].href) as src:
+                    arr[band][:,i] = np.array([x for x in src.sample(coords)]).flatten()
+            i += 1
+            prg.update()
+        prg.close()
+        api = None
+        all_arr.append(arr)
+
+    
+    days = np.array([(x-dates[0]).days for x in dates])
+    final_arr = {}
+    for b in band_list:
+        final_arr[b] = np.hstack([x[b] for x in all_arr])
+    
+    band_list.remove('SCL')
+    final_arr['SCL'] = np.invert(np.isin(final_arr['SCL'],[0,1,3,8,9,10]))
+    for i in tqdm.tqdm(range(ds.shape[0]), desc="Interpolating time series"):
+        x = days[final_arr['SCL'][i,:]]
+        for b in band_list:
+            y = final_arr[b][i,final_arr['SCL'][i,:]]
+            f = interp1d(x,y, fill_value='extrapolate')
+            final_arr[b][i,:] = f(days)
+    
+    _, ii = np.unique(days, return_index=True)
+    for b in band_list:
+        final_arr[b] = final_arr[b][:,ii]
+
+    #[os.remove(x) for x in fn]
+    os.remove(stats_fn)
+    os.remove(samples_fn)
+
+    return final_arr, dates[ii], list(ds[field_name.lower()])
 
 class S2PlanetaryTilePipeline(S2Sen2CorTilePipeline):
     NAME = 'S2-L2A-SEN2COR-PLANETARY'

VHR/vhrbase.py

@@ -228,3 +228,48 @@ class SPOT67RasterPipeline:
             self.types = self.types[:pipe_length]
             self.out_p = self.out_p[:pipe_length]
         return out
+
+class PleiadesOrthoPipeline(SPOT67RasterPipeline):
+    # BEGIN PLEIADES (ORT) VHR PROTOTYPE
+    MS_FOLDER_PATTERN = 'IMG_PHR*_MS_*'
+    PAN_FOLDER_PATTERN = 'IMG_PHR*_P_*'
+    TILE_PATTERN = '*_R*C*.TIF'
+    NDT = 0.0
+    REF_TYPE = otb.ImagePixelType_uint16
+    MS_LABEL = 'MS'
+    OPT_IN = ''
+
+class PleiadesOrthoAutoMosaicPipeline:
+    def __init__(self, mosaic_root, pattern, out_fld):
+        lst = sorted(glob.glob(os.path.join(mosaic_root, pattern)))
+        lst = [x for x in lst if os.path.isdir(x)]
+        self.scene_list = [PleiadesOrthoPipeline(x) for x in lst]
+        self.processed_scenes = []
+        self.out_fld = out_fld
+
+    def process_scenes(self, roi=None, compress=False):
+        for s in self.scene_list:
+            s.to_toa()
+            if roi is not None and os.path.exists(roi):
+                s.clip(roi)
+            s.write_outputs(self.out_fld, update_pipe=True, compress=compress)
+            s.pansharp()
+            self.processed_scenes.extend(s.write_outputs(self.out_fld, compress=compress))
+        return
+    
+    def generate_mosaic(self, vrt=True):
+        if len(self.processed_scenes) > 0:
+            lst = self.processed_scenes
+            fn = lst[0].split('_')
+            fn[-3] = 'MOSAIC'
+            fn = '_'.join(fn)
+            if vrt is True:
+                vrtopt = gdal.BuildVRTOptions()
+                fn = fn.replace('.TIF', '.vrt')
+                gdal.BuildVRT(fn, lst, options=vrtopt)
+            else:
+                print('Not implemented yet.')
+            return fn
+        else:
+            print('Preprocess scenes first.')
+            return None
\ No newline at end of file

Workflows/operations.py

@@ -32,6 +32,24 @@ def preprocess_spot67(in_fld, out_fld, dem_fld, geoid_file, skip_ps, compress,
         sp.rigid_align(align_to, this_band=align_using_band-1, ref_band=align_to_band-1)
     return sp.write_outputs(out_fld, compress=compress)
 
+def preprocess_pleiades(in_fld, out_fld, skip_ps, compress, clip, 
+                        align_to=None, align_to_band=3, align_using_band=1):
+    sp = VHR.vhrbase.PleiadesOrthoPipeline(in_fld)
+    sp.to_toa()
+    if clip is not None and os.path.exists(clip):
+        sp.clip(clip)
+    if not skip_ps:
+        sp.write_outputs(out_fld, update_pipe=True, compress=compress)
+        sp.pansharp()
+    if align_to is not None and os.path.exists(align_to):
+        sp.rigid_align(align_to, this_band=align_using_band-1, ref_band=align_to_band-1)
+    return sp.write_outputs(out_fld, compress=compress)
+
+def preprocess_pleiades_mosaic(in_fld, out_fld, pattern):
+    sp = VHR.vhrbase.PleiadesOrthoAutoMosaicPipeline(in_fld, pattern, out_fld)
+    sp.process_scenes()
+    return sp.generate_mosaic()
+
 def preprocess_s2(in_fld, out_fld, output_dates_file=None, roi=None,
                   align_to=None, align_to_band=3, align_using_band=3, provider='theia'):
     S2Processor = None

moringa.py

@@ -48,6 +48,23 @@ def main(args):
     vhrprep.add_argument("--skip_ps", help="Skip pansharpening step.", action='store_true')
     vhrprep.add_argument("--compress", help="Use lossless compression on outputs.", action='store_true')
 
+    plprep = subpar.add_parser("preprocess_pleiades", help="Perform baseline pre-processing of a Pleaides ortho scene.",
+                                formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    plprep.add_argument("fld", type=str, help="Path to the root folder containing PAN and MS subfolders.")
+    plprep.add_argument("out_fld", type=str, help="Path to the output folder for preprocessed images.")
+    plprep.add_argument("--clip", type=str, default=None, help="Path to a vector file for clipping.")
+    plprep.add_argument("--align_to", type=str, default=None, help="Path to a reference image to which the image must be aligned (rigid).")
+    plprep.add_argument("--align_to_band", type=int, default=3, help="Band of reference image used for alignment.")
+    plprep.add_argument("--align_using_band", type=int, default=1, help="Band of current image used for alignment.")
+    plprep.add_argument("--skip_ps", help="Skip pansharpening step.", action='store_true')
+    plprep.add_argument("--compress", help="Use lossless compression on outputs.", action='store_true')
+
+    plmprep = subpar.add_parser("preprocess_pleiades_mosaic", help="Perform baseline pre-processing and auto-mosaic of Pleaides ortho scenes.",
+                                formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    plmprep.add_argument("fld", type=str, help="Path to the root folder containing Pleaides acquisitions.")
+    plmprep.add_argument("pattern", type=str, help="File pattern to collect scene folders within in_fld.")
+    plmprep.add_argument("out_fld", type=str, help="Path to the output folder for preprocessed images.")
+
     s1prepr = subpar.add_parser("preprocess_s1", help="Performs Moringa preset time series preprocessing for Sentinel-1.",
                                 formatter_class=argparse.ArgumentDefaultsHelpFormatter)
     s1prepr.add_argument("in_folder", type=str, help="Path to the folder containing (de-zipped) S1 images.")
@@ -125,6 +142,13 @@ def main(args):
     if arg.cmd == "preprocess_spot67":
         preprocess_spot67(arg.fld, arg.out_fld, arg.dem_fld, arg.geoid, arg.skip_ps, arg.compress,
                           arg.clip, arg.align_to, arg.align_to_band, arg.align_using_band)
+    
+    if arg.cmd == "preprocess_pleiades":
+        preprocess_pleiades(arg.fld, arg.out_fld, arg.skip_ps, arg.compress, arg.clip, 
+                          arg.align_to, arg.align_to_band, arg.align_using_band)
+        
+    if arg.cmd == "preprocess_pleiades_mosaic":
+        preprocess_pleiades_mosaic(arg.fld, arg.out_fld, arg.pattern)
 
     if arg.cmd == "preprocess_s1":
         preprocess_s1(arg.in_folder, arg.roi, arg.out_folder, arg.dem_fld, arg.geoid, arg.direction, arg.satellite,