WIP: upscaling code to work on big areas (still under test).

Learning/ObjectBased.py

@@ -94,7 +94,7 @@ class ObjectBasedClassifier:
             }
         return models, summary, results
 
-    def classify(self, model):
+    def classify(self, model, output_file=None, compress='NONE'):
         if isinstance(model, list):
             for t, L, X in self.obia_base.tiled_data(normalize=[self.training_base['perc2'],
                                                                 self.training_base['perc98']]):
@@ -103,12 +103,12 @@ class ObjectBasedClassifier:
                     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)
+                self.obia_base.populate_map(t, L, c, output_file, compress)
         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)
+                self.obia_base.populate_map(t, L, c, output_file, compress)
         return
 
 #TEST CODE
@@ -128,6 +128,5 @@ def run_test():
     m,s,r = obc.train_RF(100)
     print(s)
     print('Performing Classification...')
-    obc.classify(m)
-    obc.obia_base.export_raster_map('/DATA/Benin/OBSYDYA_data/MORINGAv2/firstmap.tif')
+    obc.classify(m, '/DATA/Benin/OBSYDYA_data/MORINGAv2/firstmap.tif')
     return obc

OBIA/OBIABase.py

 import os.path
 import sys
 import rasterio as rio
+from rasterio.windows import Window
 from rasterio.features import shapes
 import numpy as np
 import geopandas as gpd
@@ -27,6 +28,9 @@ class OBIABase:
                  ref_data=None, ref_id_field='id', ref_class_field='class',
                  nominal_tile_size=5000, background_label=0):
         self.object_layer = object_layer
+        with rio.open(self.object_layer, 'r') as src:
+            self.crs = src.crs
+            self.transform = src.transform
         self.background_label = background_label
         if type(nominal_tile_size) is tuple:
             self.nominal_tile_size = nominal_tile_size
@@ -110,15 +114,15 @@ class OBIABase:
 
     def generate_tile_info(self, nominal_tile_size, background_label):
         in_seg = to_otb_pipeline(self.object_layer)
-        W, H = in_seg.GetImageSize('out')
+        self.W, self.H = in_seg.GetImageSize('out')
         r = otb.itkRegion()
         obj_bbox = {}
-        rw, rh = range(0, W, nominal_tile_size[0]), range(0, H, nominal_tile_size[1])
+        rw, rh = range(0, self.W, nominal_tile_size[0]), range(0, self.H, nominal_tile_size[1])
         with tqdm(total=len(rw)*len(rh), desc='Gen. tile info') as pb:
             for i in rw:
                 for j in rh:
                     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)
+                    r['size'][0], r['size'][1] = min(nominal_tile_size[0], self.W-i), min(nominal_tile_size[1], self.H-j)
                     in_seg.PropagateRequestedRegion('out', r)
                     tile = in_seg.GetImageAsNumpyArray('out')
                     rp = regionprops(tile.astype(np.uint32))
@@ -135,7 +139,7 @@ class OBIABase:
             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]]
 
-        tx, ty = int(W / nominal_tile_size[0]) + 1, int(H / nominal_tile_size[1]) + 1
+        tx, ty = int(self.W / nominal_tile_size[0]) + 1, int(self.H / nominal_tile_size[1]) + 1
 
         obj_to_tile = {}
         tiles = dict.fromkeys(range(tx * ty))
@@ -154,8 +158,6 @@ 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):
@@ -256,7 +258,7 @@ class OBIABase:
 
     def populate_ref_db(self):
         assert(self.ref_db is not None)
-        for t in tqdm(self.tiles.keys(), desc="Comp. zonal stats on DB", total=len(self.tiles)):
+        for t in tqdm(self.tiles.keys(), desc="Zonal stats on DB", total=len(self.tiles)):
             for rf, rs, rv in zip(self.raster_files, self.raster_stats, self.raster_var_names):
                 k, v = self.compute_stats_on_tile(t, rf, rs, on_ref=True)
                 if k is not None:
@@ -307,9 +309,7 @@ class OBIABase:
                 X = (X - normalize[0]) / (normalize[1] - normalize[0])
             yield tilenum,L,X
 
-    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)
+    def populate_map(self, tilenum, obj_id, classes, output_file=None, compress='NONE'):
         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)
@@ -318,8 +318,21 @@ class OBIABase:
         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]
+        if output_file is None:
+            if self.output_map is None:
+                self.output_map = np.zeros((self.H,self.W), dtype=np.int)
+            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]
+        else:
+            if tilenum == 0 and os.path.exists(output_file):
+                os.remove(output_file)
+            mode = 'w+' if not os.path.exists(output_file) else 'r+'
+            with rio.open(output_file, mode, driver='GTiff', width=self.W, height=self.H, crs=self.crs,
+                          transform=self.transform, count=1, dtype=np.int16, compress=compress) as dst:
+                win = Window(self.tiles[tilenum][0], self.tiles[tilenum][1],
+                             self.tiles[tilenum][2], self.tiles[tilenum][3])
+                x = dst.read(1, window=win)
+                dst.write(x + tmp[tile_obj], window=win, indexes=1)
         return
 
     def export_raster_map(self, output_file, compress=False):
@@ -336,6 +349,10 @@ class OBIABase:
         obj.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint16)
         obj.ExecuteAndWriteOutput()
 
+    def populate_export_raster_map(self, obj_id, classes, output_file, compress='NONE'):
+        for tn in tqdm(self.tiles.keys(), desc="Writing output map"):
+            self.populate_map(tn,obj_id,classes,output_file,compress)
+
     def true_pred_bypixel(self, labels, predicted_classes):
         pred_c = np.zeros(np.max(self.ref_db['orig_label']).astype(int)+1)
         pred_c[labels] = predicted_classes