WIP: upscaling code to work on big areas.

Common/otb_numpy_proc.py

@@ -194,7 +194,7 @@ def check_function_padding(function, test_tile_size=(100,100)):
 
 def to_otb_pipeline(obj):
     if isinstance(obj, str):
-        ppl = otb.Registry.CreateApplication('ExtractROI')
+        ppl = otb.Registry.CreateApplicationWithoutLogger('ExtractROI')
         ppl.SetParameterString('in',obj)
         ppl.Execute()
     elif isinstance(obj,otb.Application):

Learning/ObjectBased.py

@@ -94,13 +94,41 @@ class ObjectBasedClassifier:
             }
         return models, summary, results
 
+    def classify(self, model):
+        if isinstance(model, list):
+            for t, L, X in self.obia_base.tiled_data(normalize=[self.training_base['perc2'],
+                                                                self.training_base['perc98']]):
+                prob = []
+                for m in model:
+                    prob.append(m.predict_proba(X))
+                prob = np.prod(prob, axis=0)
+                c = model[0].classes_[np.argmax(prob, axis=1)]
+                self.obia_base.populate_map(t, L, c)
+        else:
+            for t,L,X in self.obia_base.tiled_data(normalize=[self.training_base['perc2'],
+                                                              self.training_base['perc98']]):
+                c = model.predict(X)
+                self.obia_base.populate_map(t, L, c)
+        return
+
 #TEST CODE
 def run_test():
+    '''
     obc = ObjectBasedClassifier('/DATA/Moringa_Sample/Parakou/output/segmentation/segmentation.tif',
                                 '/DATA/Moringa_Sample/Parakou/input/REF/ref.shp',
                                 ['/DATA/Moringa_Sample/Parakou/output/S2_processed/T31PDL/*/*FEAT.tif'],
                                 ['/DATA/Moringa_Sample/Parakou/input/THR/THR_SPOT6.tif'])
+    '''
+    print('Initializing the OBIA training base...')
+    obc = ObjectBasedClassifier('/DATA/Benin/OBSYDYA_data/MORINGA/SEGMENTATION/segmentation.tif',
+                                '/DATA/Benin/OBSYDYA_data/MORINGA/reference/BD_OBSYDYA_2022_ParakouNdali_v0.2.shp',
+                                ['/DATA/Benin/OBSYDYA_data/MORINGA/basefolder/FEAT/S2_THEIA_FEAT/S2_THEIA_MOSAIC_*.tif'],
+                                glob.glob('/DATA/Benin/OBSYDYA_data/MORINGA/ext_features'))
     obc.gen_k_folds(5)
-    m,s,r = obc.train_RF(100)
+    print('Performing Training and Cross-Validation...')
+    m,s,r = obc.train_RF(400)
     print(s)
-    return
+    print('Performing Classification...')
+    obc.classify(m)
+    obc.obia_base.export_raster_map('/DATA/Benin/OBSYDYA_data/MORINGAv2/firstmap.tif')
+    return obc

OBIA/OBIABase.py

@@ -48,7 +48,11 @@ class OBIABase:
         self.raster_groups = []
         self.n_vars = 0
 
+        # Output raster
+        self.output_map = None
+
     def init_ref_db(self, vector_file, id_field, class_field):
+        print("Enter init")
         ras_id = otb.Registry.CreateApplication('Rasterization')
         ras_id.SetParameterString('in', vector_file)
         ras_id.SetParameterString('im', self.object_layer)
@@ -56,6 +60,7 @@ class OBIABase:
         ras_id.SetParameterString('mode.attribute.field', id_field)
         ras_id.Execute()
         ids = ras_id.GetImageAsNumpyArray('out')
+        print("First raster")
 
         ras_cl = otb.Registry.CreateApplication('Rasterization')
         ras_cl.SetParameterString('in', vector_file)
@@ -64,15 +69,18 @@ class OBIABase:
         ras_cl.SetParameterString('mode.attribute.field', class_field)
         ras_cl.Execute()
         cls = ras_cl.GetImageAsNumpyArray('out')
+        print("Second raster")
 
         in_seg = to_otb_pipeline(self.object_layer)
         obj = in_seg.GetImageAsNumpyArray('out')
+        print("Load segmentation")
 
         mask = ids>0
         # This was to relabel object being disconnected by intersecting the GT
         # Re-enable if needed...
         # self.ref_obj_layer = label(obj * mask, connectivity=1).astype(np.uint32)
         self.ref_obj_layer = (obj * mask).astype(np.uint32)
+        print("Before Regionprops")
         rp = regionprops(self.ref_obj_layer, intensity_image=np.dstack((obj,ids,cls)).astype(np.int))
         self.ref_db = pd.DataFrame(data=[np.insert(o.intensity_min,0,o.area) for o in rp],
                                    columns=['area','orig_label','polygon_id','class'],
@@ -80,11 +88,35 @@ class OBIABase:
         return
 
     def generate_tile_info(self, nominal_tile_size, background_label):
+        print('Enter generate_tile')
         in_seg = to_otb_pipeline(self.object_layer)
-        full = in_seg.GetImageAsNumpyArray('out')
-        rp = regionprops(full.astype(np.uint32))
-
         W, H = in_seg.GetImageSize('out')
+        r = otb.itkRegion()
+        obj_bbox = {}
+        for i in range(0, W, nominal_tile_size[0]):
+            for j in range(0, H, nominal_tile_size[1]):
+                r['index'][0], r['index'][1] = i, j
+                r['size'][0], r['size'][1] = min(nominal_tile_size[0], W-i), min(nominal_tile_size[1], H-j)
+                in_seg.PropagateRequestedRegion('out', r)
+                tile = in_seg.GetImageAsNumpyArray('out')
+                rp = regionprops(tile.astype(np.uint32))
+                for o in rp:
+                    bbox = np.array([o.bbox[0]+j, o.bbox[1]+i, o.bbox[2]+j, o.bbox[3]+i])
+                    if o.label in obj_bbox.keys():
+                        obj_bbox[o.label].append(bbox)
+                    else:
+                        obj_bbox[o.label] = [bbox]
+
+        for o in obj_bbox.keys():
+            obj_bbox[o] = np.array(obj_bbox[o])
+            mn, mx = np.min(obj_bbox[o], axis=0), np.max(obj_bbox[o], axis=0)
+            obj_bbox[o] = [mn[0], mn[1], mx[2], mx[3]]
+
+        print("Finished BBOX.")
+        print("BBOX:")
+        print(obj_bbox)
+        input("Pippo...")
+
         tx, ty = int(W / nominal_tile_size[0]) + 1, int(H / nominal_tile_size[1]) + 1
 
         obj_to_tile = {}
@@ -104,6 +136,8 @@ class OBIABase:
 
         tiles = {k: [v[0],v[1],v[2]-v[0],v[3]-v[1]] for k,v in tiles.items() if not np.isinf(v[0])}
 
+        self.W, self.H = W, H
+
         return tiles, obj_to_tile
 
     def compute_stats_on_tile(self, tile_num, input_image, stats=[OStats.MEAN], on_ref=False):
@@ -246,15 +280,41 @@ class OBIABase:
         vars = [item for sublist in self.raster_var_names for item in sublist]
         for tilenum in self.tiles.keys():
             tile_data = self.get_full_stats_on_tile(tilenum)
-            L = tile_data.index.to_numpy(dtype=np.int)
+            L = tile_data.index.to_numpy(dtype=np.int32)
             X = tile_data[vars].to_numpy()
             if normalize is not None:
                 X = (X - normalize[0]) / (normalize[1] - normalize[0])
-            yield L,X
+            yield tilenum,L,X
 
-    def populate_map(self, tilenum, labels):
+    def populate_map(self, tilenum, obj_id, classes):
+        if self.output_map is None:
+            self.output_map = np.zeros((self.H,self.W), dtype=np.int)
+        r = otb.itkRegion()
+        r['index'][0], r['index'][1], r['size'][0], r['size'][1] = self.tiles[tilenum]
+        obj = to_otb_pipeline(self.object_layer)
+        obj.PropagateRequestedRegion('out', r)
+        clip_obj = obj.ExportImage('out')
+        tile_obj = np.squeeze(clip_obj['array']).astype(np.uint32)
+        tmp = np.zeros(np.max(tile_obj) + 1, dtype=int)
+        tmp[obj_id] = classes
+        self.output_map[self.tiles[tilenum][1]:self.tiles[tilenum][1]+self.tiles[tilenum][3],
+        self.tiles[tilenum][0]:self.tiles[tilenum][0]+self.tiles[tilenum][2]] += tmp[tile_obj]
         return
 
+    def export_raster_map(self, output_file, compress=False):
+        assert (self.output_map is not None)
+        if not os.path.exists(os.path.dirname(output_file)):
+            os.makedirs(os.path.dirname(output_file))
+        obj = to_otb_pipeline(self.object_layer)
+        arr = obj.ExportImage('out')
+        arr['array'] = self.output_map.astype(np.float32)
+        obj.ImportImage('in', arr)
+        if compress:
+            output_file = output_file + "?gdal:co:compress=deflate"
+        obj.SetParameterString('out', output_file)
+        obj.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint16)
+        obj.ExecuteAndWriteOutput()
+
     def true_pred_bypixel(self, labels, predicted_classes):
         pred_c = np.zeros(np.max(self.ref_obj_layer)+1)
         pred_c[labels] = predicted_classes

moringa.py

@@ -91,7 +91,7 @@ def main(args):
     vhrprep.add_argument("--align_to_band", type=int, default=3, help="Band of reference image used for alignment.")
     vhrprep.add_argument("--align_using_band", type=int, default=1, help="Band of current image used for alignment.")
     vhrprep.add_argument("--skip_ps", help="Skip pansharpening step.", action='store_true')
-    vhrprep.add_argument("--compress", help="Use lossless compresion on outputs.", action='store_true')
+    vhrprep.add_argument("--compress", help="Use lossless compression on outputs.", action='store_true')