OBIA/OBIABase.py

@@ -377,14 +377,14 @@ class OBIABase:
             self.populate_map(tn,obj_id,classes,output_file,compress)
 
     def true_pred_bypixel(self, labels, predicted_classes, class_field='class'):
-        pred_c = np.zeros(np.max(self.ref_db['orig_label']).astype(int)+1)
+        pred_c = np.zeros(int(np.max(self.ref_db['orig_label']))+1)
         pred_c[labels] = predicted_classes
         support = []
         for tn, t in self.tiled_objects(on_ref=True):
             support.append(t[np.isin(t, labels)])
         support = np.concatenate(support)
         pred = pred_c[support]
-        true_c = np.zeros(np.max(self.ref_db['orig_label']).astype(int)+1)
+        true_c = np.zeros(int(np.max(self.ref_db['orig_label']))+1)
         # ATTENTION: works if "labels" is sorted (as provided by get_reference_...)
         true_c[labels] = self.ref_db.loc[self.ref_db['orig_label'].isin(labels),class_field].to_numpy(dtype=int)
         true = true_c[support]

TimeSeries/landsat_planetary.py

+from TimeSeries.s2theia import *
 from Common.geotools import get_query_bbox
 from pystac_client import Client
 import planetary_computer as PC
@@ -52,4 +53,133 @@ def fetch(shp, dt, out_fld, auth=None, band_list=None):
             with rasterio.open(ofn, "w", **out_meta) as dest:
                 dest.write(out_img)
             prg.update()
-    prg.close()
\ No newline at end of file
+    prg.close()
+
+
+class Landsat89PlanetaryTilePipeline(S2TheiaTilePipeline):
+    # --- BEGIN SENSOR PROTOTYPE ---
+
+    NAME = 'LANDSAT89-PLANETARY'
+    REF_TYPE = otb.ImagePixelType_uint16
+    MSK_TYPE = otb.ImagePixelType_uint8
+    PTRN_dir = 'LC0*_L2SP*'
+    PTRN_ref = '_B'
+    B2_name = 'B1'
+    PTRN_30m = ['*_SR_B1.TIF', '*_SR_B2.TIF', '*_SR_B3.TIF',
+                '*_SR_B4.TIF', '*_SR_B5.TIF', '*_SR_B6.TIF', 
+                '*_SR_B7.TIF']
+    PTRN_msk = ['*_QA_PIXEL.TIF']
+    PTRN_ful = PTRN_30m
+    MASK_exp = 'im1b1!=21824 && im1b1!=21888 && im1b1!=30048'
+    FEAT_exp = {
+        'B1': 'im1b1',
+        'B2': 'im1b2',
+        'B3': 'im1b3',
+        'B4': 'im1b4',
+        'B5': 'im1b5',
+        'B6': 'im1b6',
+        'B7': 'im1b7',
+        'NDVI': '(im1b5-im1b4)/(im1b5+im1b4+1e-6)',
+        'NDWI': '(im1b3-im1b5)/(im1b3+im1b5+1e-6)',
+        'BRI': 'sqrt(' + '+'.join(['im1b%d*im1b%d' % (i, i) for i in range(1, 8)]) + ')',
+        'MNDWI': '(im1b3-im1b7)/(im1b3+im1b7+1e-6)',
+        'SWNDVI': '(im1b7-im1b5)/(im1b7+im1b5+1e-6)',
+    }
+    NDT = 0
+
+    @classmethod
+    def _check(cls,x):
+        return os.path.isdir(x) and os.path.basename(x).startswith(cls.PTRN_dir.split('*')[0]) and 'L2SP' in os.path.basename(x)
+
+    @classmethod
+    def _img_id(cls,x):
+        if cls._check(x):
+            return os.path.basename(x)
+        else:
+            return None
+
+    @classmethod
+    def _img_date(cls,x):
+        if cls._check(x):
+            return os.path.basename(x).split('_')[3]
+        else:
+            return None
+
+    @classmethod
+    def _tile_id(cls,x):
+        if cls._check(x):
+            return os.path.basename(x).split('_')[2]
+        else:
+            return None
+
+    def _process_mask(self, msks):
+        msk_pipe = otb.Registry.CreateApplication('BandMath')
+        msk_pipe.AddImageToParameterInputImageList('il', msks[0].GetParameterOutputImage('out'))
+        msk_pipe.SetParameterString('exp', self.MASK_exp)
+        msk_pipe.Execute()
+        return msk_pipe
+    
+    # ---- END SENSOR PROTOTYPE ----
+    
+    def preprocess(self):
+        proc_idx = self.out_idx.copy()
+        self.out_idx = []
+        num_of_band_inputs = len(self.PTRN_ful)
+        num_of_mask_inputs = len(self.PTRN_msk)
+        for t in proc_idx[::num_of_band_inputs+num_of_mask_inputs]:
+            idx_to_stack = []
+            cct = otb.Registry.CreateApplication('ConcatenateImages')
+            [cct.AddImageToParameterInputImageList('il', self.pipe[j].GetParameterOutputImage('out')) for j in range(t, t+num_of_band_inputs)]
+            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+num_of_band_inputs:t+num_of_band_inputs+num_of_mask_inputs]
+            msk_pipe = self._process_mask(msks)
+            fn = self.files[t].replace(self.B2_name, 'BINARY_MASK')
+            ty = self.MSK_TYPE
+            self.append(msk_pipe, fn, ty, 'out', is_output=True)
+
+class Landsat7PlanetaryTilePipeline(Landsat89PlanetaryTilePipeline):
+    NAME = 'LANDSAT7-PLANETARY'
+    PTRN_dir = 'LE07_L2SP*'
+    PTRN_30m = ['*_SR_B1.TIF', '*_SR_B2.TIF', '*_SR_B3.TIF',
+                '*_SR_B4.TIF', '*_SR_B5.TIF', '*_SR_B7.TIF']
+    PTRN_ful = PTRN_30m
+    MASK_exp = 'im1b1!=5440 && im1b1!=5504 && im1b1!=13664'
+    FEAT_exp = {
+        'B1': 'im1b1',
+        'B2': 'im1b2',
+        'B3': 'im1b3',
+        'B4': 'im1b4',
+        'B5': 'im1b5',
+        'B7': 'im1b6',
+        'NDVI': '(im1b4-im1b3)/(im1b4+im1b3+1e-6)',
+        'NDWI': '(im1b2-im1b4)/(im1b2+im1b4+1e-6)',
+        'BRI': 'sqrt(' + '+'.join(['im1b%d*im1b%d' % (i, i) for i in range(1, 7)]) + ')',
+        'MNDWI': '(im1b2-im1b6)/(im1b2+im1b6+1e-6)',
+        'SWNDVI': '(im1b6-im1b4)/(im1b6+im1b4+1e-6)',
+    }
+
+class Landsat45PlanetaryTilePipeline(Landsat7PlanetaryTilePipeline):
+    NAME = 'LANDSAT45-PLANETARY'
+    PTRN_dir = 'LT0*_L2SP*'
+    
+class Landsat89PlanetaryPipeline(S2TheiaPipeline):
+
+    S2TilePipeline = Landsat89PlanetaryTilePipeline
+    _check = S2TilePipeline._check
+    _tile_id = S2TilePipeline._tile_id
+
+class Landsat7PlanetaryPipeline(S2TheiaPipeline):
+
+    S2TilePipeline = Landsat7PlanetaryTilePipeline
+    _check = S2TilePipeline._check
+    _tile_id = S2TilePipeline._tile_id
+
+class Landsat45PlanetaryPipeline(S2TheiaPipeline):
+
+    S2TilePipeline = Landsat45PlanetaryTilePipeline
+    _check = S2TilePipeline._check
+    _tile_id = S2TilePipeline._tile_id

TimeSeries/s2theia.py

@@ -163,9 +163,9 @@ class S2TheiaTilePipeline:
     def _check_roi(cls, x, roi, min_surf=0.0, temp_fld='/tmp'):
         if cls._check(x):            
             if os.path.isdir(x):
-                bnd = glob.glob(os.path.join(x, cls.PTRN_20m[0]))[0]
+                bnd = glob.glob(os.path.join(x, cls.PTRN_ful[0]))[0]
             elif os.path.splitext(x)[-1] == '.zip':
-                idf = cls.PTRN_20m[0].replace('*', '')
+                idf = cls.PTRN_ful[0].replace('*', '')
                 arch = zipfile.ZipFile(x)
                 fid = [name for name in arch.namelist() if idf in name]
                 fle = arch.read(fid[0])
@@ -340,7 +340,7 @@ class S2TheiaTilePipeline:
             i += 1
         to_merge.append(curr)
 
-        T = 10 + len(self.PTRN_msk)
+        T = len(self.PTRN_ful) + len(self.PTRN_msk)
         new_dates = []
         proc_idx = self.out_idx.copy()
         self.out_idx = []
@@ -351,18 +351,18 @@ class S2TheiaTilePipeline:
                     idx = T*l[0]+k
                     self.append(self.pipe[idx], self.files[idx], self.types[idx], 'out', is_output=True)
             else:
-                for k in range(10):
+                for k in range(len(self.PTRN_ful)):
                     mos = otb.Registry.CreateApplication('Mosaic')
                     for i in l:
                         mos.AddImageToParameterInputImageList('il', self.pipe[T*i+k].GetParameterOutputImage(self.out_p[T*i+k]))
                     mos.SetParameterInt('nodata', self.NDT)
                     mos.Execute()
                     self.append(mos, self.files[T*i+k], self.types[T*i+k], 'out', is_output=True)
-                for k in range(10,T):
+                for k in range(len(self.PTRN_ful),T):
                     bm = otb.Registry.CreateApplication('BandMath')
                     for i in l:
                         bm.AddImageToParameterInputImageList('il', self.pipe[T*i+k].GetParameterOutputImage(self.out_p[T*i+k]))
-                    bm.SetParameterString('exp', self.MERG_msk[k-10] + '({})'.format(','.join(['im{}b1'.format(w+1) for w in range(len(l))])))
+                    bm.SetParameterString('exp', self.MERG_msk[k-len(self.PTRN_ful)] + '({})'.format(','.join(['im{}b1'.format(w+1) for w in range(len(l))])))
                     bm.Execute()
                     self.append(bm, self.files[T * i + k], self.types[T * i + k], 'out', is_output=True)
             new_dates.append(self.input_dates[l[-1]])
@@ -385,7 +385,7 @@ class S2TheiaTilePipeline:
         return cc
 
     def set_input_epsg(self):
-        f, _ = self.get_file(self.image_list[0], self.PTRN_20m[0])
+        f, _ = self.get_file(self.image_list[0], self.PTRN_ful[0])
         ds = gdal.Open(f)
         self.input_epsg = osr.SpatialReference(wkt=ds.GetProjection()).GetAuthorityCode('PROJCS')
         ds = None
@@ -417,10 +417,10 @@ class S2TheiaTilePipeline:
         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]:
+        for t in proc_idx[::len(self.PTRN_ful)+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]:
+            for k in range(t,t+len(self.PTRN_ful)):
+                if k % (len(self.PTRN_ful) + 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'))
@@ -442,7 +442,7 @@ class S2TheiaTilePipeline:
             ty = self.REF_TYPE
             self.append(cct, fn, ty, 'out', is_output=True)
 
-            msks = self.pipe[t+10:t+10+num_of_mask_inputs]
+            msks = self.pipe[t+len(self.PTRN_ful):t+len(self.PTRN_ful)+num_of_mask_inputs]
             msk_pipe = self._process_mask(msks)
             for app in msk_pipe[:-1]:
                 self.append(app)
@@ -554,8 +554,8 @@ class S2TheiaTilePipeline:
         proc_idx = self.out_idx.copy()
         self.out_idx = []
         img = [self.pipe[t] for t in range(self.pipe_start+match_band,
-                                           self.pipe_start+(10+len(self.PTRN_msk))*len(self.input_dates),
-                                           10+len(self.PTRN_msk))]
+                                           self.pipe_start+(len(self.PTRN_ful)+len(self.PTRN_msk))*len(self.input_dates),
+                                           len(self.PTRN_ful)+len(self.PTRN_msk))]
         msk = [self.pipe[t] for t in proc_idx[1::2]]
         RX, RY = img[0].GetImageSpacing(out_param)
         ref_idx = None
@@ -662,7 +662,7 @@ class S2TheiaTilePipeline:
         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.SetParameterInt('comp', len(self.PTRN_ful))
         gf.SetParameterString('it', 'linear')
         gf.SetParameterString('id', os.path.join(self.folder, self.tile_id + '_indates.txt'))
         gf.SetParameterString('od', output_dates)
@@ -670,7 +670,7 @@ class S2TheiaTilePipeline:
             gf.Execute()
             gf_out = gf
         else:
-            gf_fn = self.temp_fld + os.sep + 'SENTINEL2_' + self.tile_id + '_GAPFILLED_FULL.tif'
+            gf_fn = self.temp_fld + os.sep + '{}_'.format(self.NAME) + self.tile_id + '_GAPFILLED_FULL.tif'
             if not os.path.exists(gf_fn):
                 gf.SetParameterString("out", gf_fn)
                 gf.ExecuteAndWriteOutput()
@@ -679,14 +679,14 @@ class S2TheiaTilePipeline:
 
         t = 1
         for d in od:
-            ch_list = ['Channel%d' % i for i in range(t,t+10)]
-            t += 10
+            ch_list = ['Channel%d' % i for i in range(t,t+len(self.PTRN_ful))]
+            t += len(self.PTRN_ful)
             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
+            dn = '{}_'.format(self.NAME) + 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)
@@ -726,7 +726,7 @@ class S2TheiaTilePipeline:
                 i += 1
             bm.SetParameterString('exp', '{' + ';'.join(expr) + '}')
             bm.Execute()
-            fn = self.tile_id + os.sep + 'SENTINEL2_' + self.tile_id + '_GAPFILL_' + feat + '.tif'
+            fn = self.tile_id + os.sep + '{}_'.format(self.NAME) + self.tile_id + '_GAPFILL_' + feat + '.tif'
             ty = otb.ImagePixelType_float
             self.append(bm, fn, ty, 'out', is_output=True)
 
@@ -892,7 +892,7 @@ class S2TheiaPipeline:
                 out_mos = []
                 vrtopt = gdal.BuildVRTOptions()
                 for i in range(n_dates):
-                    fn = out_fld + os.sep + 'SENTINEL2_MOSAIC_GAPFILL_' + self.tiles[0].output_dates[i] + '_' + stack_name + '.vrt'
+                    fn = out_fld + os.sep + '{}_MOSAIC_GAPFILL_'.format(self.tiles[0].NAME) + 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)
@@ -916,7 +916,7 @@ class S2TheiaPipeline:
                 out_mos = []
                 vrtopt = gdal.BuildVRTOptions()
                 for i in range(len(feat_list)):
-                    fn = out_fld + os.sep + 'SENTINEL2_MOSAIC_GAPFILL_' + feat_list[i] + '.vrt'
+                    fn = out_fld + os.sep + '{}_MOSAIC_GAPFILL_'.format(self.tiles[0].NAME) + feat_list[i] + '.vrt'
                     to_mosaic = [x[i] for x in out]
                     gdal.BuildVRT(fn, to_mosaic, options=vrtopt)
                     out_mos.append(fn)

Workflows/basic.py

@@ -51,6 +51,16 @@ def process_timeseries(oroot, d, ts_lst_pkl):
             os.makedirs(ots, exist_ok=True)
             ts_lst.append(preprocess_planet(ts['path'],
                                             ots))
+        elif ts['type'].startswith('landsat'):
+            print('[MORINGA-INFO] : Preprocessing {}'.format(ts['type']))
+            ots = os.path.join(oroot, 'timeseries/' + ts['type'])
+            os.makedirs(ots, exist_ok=True)
+            ts_lst.append(preprocess_landsat(ts['path'],
+                                             ots,
+                                             roi=d['roi'],
+                                             output_dates_file=ts['output_dates_file'],
+                                             which=ts['type']))
+
         else:
             raise ValueError('TimeSeries type not yet supported.')
     with open(ts_lst_pkl, 'wb') as ts_save:
@@ -84,7 +94,9 @@ def train_valid_workflow(seg, ts_lst_pkl, d, m_file):
                                 reference_data=d['ref_db']['path'],
                                 ref_class_field=d['ref_db']['fields'])
 
-    obc.gen_k_folds(5, class_field=d['ref_db']['fields'][-1],augment=d['training']['augment_if_missing'])
+    obc.gen_k_folds(5, class_field=d['ref_db']['fields'][-1],
+                    augment=d['training']['augment_if_missing'], 
+                    min_samples_per_class=d['training']['augment_if_missing'])
 
     if 'export_training_base' in d['training'].keys() and d['training']['export_training_base'] is True:
         obc.save_training_base('{}/_side/training_base.pkl'.format(os.path.join(d['output_path'], d['chain_name'])))

Workflows/operations.py

@@ -55,6 +55,25 @@ def preprocess_s2(in_fld, out_fld, output_dates_file=None, roi=None,
     return s2.extract_feature_set(out_fld, store_gapfill=True, mosaicking='vrt', align=align,
                            align_to=align_to, align_to_band=align_to_band, align_using_band=align_using_band)
 
+def preprocess_landsat(in_fld, out_fld, output_dates_file=None, roi=None, which='landsat89'):
+    LandsatProcessor = None
+    if which == 'landsat89':
+        LandsatProcessor = landsat_planetary.Landsat89PlanetaryPipeline
+    elif which == 'landsat7':
+        LandsatProcessor = landsat_planetary.Landsat7PlanetaryPipeline
+    elif which == 'landsat45':
+        LandsatProcessor = landsat_planetary.Landsat45PlanetaryPipeline
+    else:
+        raise ValueError("Unknown Landsat mission.")
+    
+    ls = LandsatProcessor(in_fld, temp_fld=out_fld, roi=roi)
+    if output_dates_file is not None:
+        ls.set_output_dates_by_file(output_dates_file)
+    else:
+        raise ValueError("Please provide path to a text file containing output dates.")
+    
+    return ls.extract_feature_set(out_fld, store_gapfill=True, mosaicking='vrt')
+
 def preprocess_planet(in_fld, out_fld):
     pl = planet_mosaics.PlanetMosaicPipeline(in_fld)
     pl.compute_features()