diff --git a/Learning/ObjectBased.py b/Learning/ObjectBased.py index 1fa7787723341ec8e7e831372e958d22c568cad4..aa7aa9636c2bdbe71b5364b66050ac91d2a8ea59 100644 --- a/Learning/ObjectBased.py +++ b/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 diff --git a/OBIA/OBIABase.py b/OBIA/OBIABase.py index 3e264e49b4b97bd6b7c715037f5a1c9e434445ef..a5e33f726873bf3fb35d1c24d8295a44694b5efe 100644 --- a/OBIA/OBIABase.py +++ b/OBIA/OBIABase.py @@ -1,6 +1,7 @@ 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