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 import pandas as pd from enum import Enum from Common.otb_numpy_proc import to_otb_pipeline from skimage.measure import regionprops, label import otbApplication as otb from tqdm import tqdm from sklearn.impute import SimpleImputer class OStats(Enum): MEAN = 'mean' STD = 'std' MIN = 'min' MAX = 'max' def intensity_std(region, intensities): # note the ddof arg to get the sample var if you so desire! return np.std(intensities[region]) class OBIABase: def __init__(self, object_layer, 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 elif type(nominal_tile_size) is int: self.nominal_tile_size = (nominal_tile_size, nominal_tile_size) else: sys.exit('Invalid nominal tile size') self.tiles, self.obj_to_tile = self.generate_tile_info(self.nominal_tile_size, self.background_label) # Init database self.ref_db = None self.ref_obj_layer_pipe = None if ref_data is not None: assert(ref_id_field is not None and ref_class_field is not None) self.init_ref_db(ref_data, ref_id_field, ref_class_field) # Computation rasters self.raster_files = [] self.raster_stats = [] self.raster_var_names = [] self.raster_groups = [] self.n_vars = 0 # Output raster self.output_maps = [] def init_ref_db(self, vector_file, id_field, class_field): if isinstance(class_field, str): class_field = [class_field] ras_id = otb.Registry.CreateApplication('Rasterization') ras_id.SetParameterString('in', vector_file) ras_id.SetParameterString('im', self.object_layer) ras_id.SetParameterString('mode', 'attribute') ras_id.SetParameterString('mode.attribute.field', id_field) ras_id.Execute() #ids = ras_id.GetImageAsNumpyArray('out') ras_cl = [] for cf in class_field: ras_cl.append(otb.Registry.CreateApplication('Rasterization')) ras_cl[-1].SetParameterString('in', vector_file) ras_cl[-1].SetParameterString('im', self.object_layer) ras_cl[-1].SetParameterString('mode', 'attribute') ras_cl[-1].SetParameterString('mode.attribute.field', cf) ras_cl[-1].Execute() in_seg = to_otb_pipeline(self.object_layer) intensity_img = otb.Registry.CreateApplication('ConcatenateImages') intensity_img.AddImageToParameterInputImageList('il', in_seg.GetParameterOutputImage('out')) intensity_img.AddImageToParameterInputImageList('il', ras_id.GetParameterOutputImage('out')) [intensity_img.AddImageToParameterInputImageList('il', rcl.GetParameterOutputImage('out')) for rcl in ras_cl] intensity_img.Execute() ref_ol = otb.Registry.CreateApplication('BandMath') ref_ol.AddImageToParameterInputImageList('il', in_seg.GetParameterOutputImage('out')) ref_ol.AddImageToParameterInputImageList('il', ras_id.GetParameterOutputImage('out')) ref_ol.SetParameterString('exp', 'im1b1 * (im2b1 > 0)') ref_ol.Execute() self.ref_obj_layer_pipe = [in_seg, ras_id, ref_ol] self.ref_db = None ''' self.ref_db = pd.DataFrame(data=[], columns=['area', 'orig_label', 'polygon_id'] + class_field, index=[]) ''' r = otb.itkRegion() for tn, t in tqdm(self.tiles.items(), desc='Init. Ref. DB', total=len(self.tiles)): r['index'][0], r['index'][1] = t[0], t[1] r['size'][0], r['size'][1] = t[2], t[3] ref_ol.PropagateRequestedRegion('out', r) intensity_img.PropagateRequestedRegion('out', r) tile_ref_ol = ref_ol.GetImageAsNumpyArray('out').astype(np.int32) tile_int_img = intensity_img.GetVectorImageAsNumpyArray('out').astype(int) rp = regionprops(tile_ref_ol, intensity_image=tile_int_img) if self.ref_db is None: self.ref_db = pd.DataFrame(data=[np.insert(o.intensity_min, 0, o.area) for o in rp if self.obj_to_tile[o.label] == tn], columns=['area', 'orig_label', 'polygon_id'] + class_field, index=[o.label for o in rp if self.obj_to_tile[o.label] == tn]) else: self.ref_db = pd.concat([ self.ref_db, pd.DataFrame( data=[np.insert(o.intensity_min, 0, o.area) for o in rp if self.obj_to_tile[o.label] == tn], columns=self.ref_db.columns, index=[o.label for o in rp if self.obj_to_tile[o.label] == tn] )] ) return def generate_tile_info(self, nominal_tile_size, background_label): in_seg = to_otb_pipeline(self.object_layer) self.W, self.H = in_seg.GetImageSize('out') r = otb.itkRegion() obj_bbox = {} 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], 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)) 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] pb.update(1) 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]] 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)) for i in range(tx * ty): tiles[i] = [np.inf, np.inf, 0, 0] for o,bb in obj_bbox.items(): if o != background_label: ix, iy = int(bb[1] / nominal_tile_size[0]), int(bb[0] / nominal_tile_size[1]) idx = ix * ty + iy obj_to_tile[o] = idx tiles[idx][0] = int(min(bb[1], tiles[idx][0])) tiles[idx][1] = int(min(bb[0], tiles[idx][1])) tiles[idx][2] = int(max(bb[3], tiles[idx][2])) tiles[idx][3] = int(max(bb[2], tiles[idx][3])) 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])} return tiles, obj_to_tile def compute_stats_on_tile(self, tile_num, input_image, stats=[OStats.MEAN], on_ref=False): r = otb.itkRegion() r['index'][0], r['index'][1] = self.tiles[tile_num][0], self.tiles[tile_num][1] r['size'][0], r['size'][1] = self.tiles[tile_num][2], self.tiles[tile_num][3] if not on_ref: 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) else: assert (self.ref_obj_layer_pipe is not None) tmp_er = otb.Registry.CreateApplication('ExtractROI') tmp_er.SetParameterInputImage('in', self.ref_obj_layer_pipe[-1].GetParameterOutputImage('out')) tmp_er.SetParameterInt('startx', r['index'][0]) tmp_er.SetParameterInt('starty', r['index'][1]) tmp_er.SetParameterInt('sizex', r['size'][0]) tmp_er.SetParameterInt('sizey', r['size'][1]) tmp_er.Execute() tile_obj = tmp_er.GetImageAsNumpyArray('out').astype(np.uint32) si = otb.Registry.CreateApplication('Superimpose') si.SetParameterString('inm', input_image) si.SetParameterString('inr', self.object_layer) si.SetParameterString('interpolator', 'nn') si.Execute() si.PropagateRequestedRegion('out', r) clip_img = si.ExportImage('out') tile_img = np.squeeze(clip_img['array']) if OStats.STD in stats: op = regionprops(tile_obj, intensity_image=tile_img, extra_properties=[intensity_std]) else: op = regionprops(tile_obj, intensity_image=tile_img) out_stats = {} for o in op: if (not on_ref and self.obj_to_tile[o.label] == tile_num) or \ (on_ref and self.obj_to_tile[self.ref_db["orig_label"][o.label]] == tile_num): out_stats[o.label] = [] if OStats.MEAN in stats: out_stats[o.label].append(o.intensity_mean) if OStats.STD in stats: out_stats[o.label].append(o.intensity_std) if OStats.MIN in stats: out_stats[o.label].append(o.intensity_min) if OStats.MAX in stats: out_stats[o.label].append(o.intensity_max) out_stats[o.label] = np.array(out_stats[o.label]).flatten(order='F') if len(out_stats) > 0: return out_stats.keys(), np.vstack(list(out_stats.values())) else: return None, None def add_raster_for_stats(self, raster_file, raster_name=None, stats=[OStats.MEAN], as_group=True): self.raster_files.append(raster_file) self.raster_stats.append(stats) if raster_name is None: raster_name = os.path.splitext(os.path.basename(raster_file))[0] with rio.open(raster_file) as ds: n_bands = ds.count ds = None to_append = [raster_name + '_{}_{}'.format(i+1,j.value) for i in range(n_bands) for j in stats] if as_group: self.raster_groups.extend(range(self.n_vars, self.n_vars+len(to_append))) self.raster_var_names.append(to_append) self.n_vars += len(to_append) return def add_raster_time_series_for_stats(self, raster_list, ts_name=None, stats=[OStats.MEAN]): #Check n_bands = [] for r in raster_list: with rio.open(r) as ds: n_bands.append(ds.count) assert all(x == n_bands[0] for x in n_bands) n_bands = n_bands[0] pos = self.n_vars for i,r in enumerate(raster_list): if ts_name is not None: raster_name = '{}_D{}'.format(ts_name,i) self.add_raster_for_stats(r, stats=stats, raster_name=ts_name, as_group=False) for i in range(n_bands): for s in range(len(stats)): self.raster_groups.append(list(range(len(stats)*i+s+pos, len(stats)*i+s+pos+len(stats)*n_bands*len(raster_list), len(stats)*n_bands))) return def get_full_stats_on_tile(self, tilenum): out_db = pd.DataFrame(index=[o for o,t in self.obj_to_tile.items() if t == tilenum]) for rf,rs,rv in zip(self.raster_files, self.raster_stats, self.raster_var_names): k,v = self.compute_stats_on_tile(tilenum, rf, rs) if k is not None: out_db.loc[k,rv] = v return out_db def populate_ref_db(self): assert(self.ref_db is not None) 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: self.ref_db.loc[k,rv] = v return def export_reference_db(self, out_vector): #TO BE CHANGED with rio.open(self.object_layer) as ds: feats = [{'properties': {'sample_id': v}, 'geometry': s} for s, v in shapes(self.ref_obj_layer.astype(np.int32), mask=self.ref_obj_layer > 0, transform=ds.transform)] crs = ds.crs.to_epsg() vds = gpd.GeoDataFrame.from_features(feats) vds = vds.set_crs(epsg=crs) vds = vds.join(self.ref_db, on="sample_id") vds.to_file(out_vector) return def get_vars(self): return [item for sublist in self.raster_var_names for item in sublist] def get_reference_db_as_training_base(self, class_field='class'): if isinstance(class_field, str): class_field = [class_field] assert(self.ref_db is not None and len(self.raster_var_names)>0) vars = self.get_vars() out = {} out['vars'] = vars out['obj_id'] = self.ref_db['orig_label'].to_numpy(dtype=int) imp = SimpleImputer(missing_values=np.nan, strategy='mean') out['X'] = imp.fit_transform(self.ref_db[vars].to_numpy()) # compute percentiles and normalize out['perc2'] = np.zeros(out['X'].shape[1]) out['perc98'] = np.zeros(out['X'].shape[1]) for g in self.raster_groups: tmp = out['X'][:,g] m,M = np.percentile(tmp, [2, 98]) if isinstance(g, list): for x in g: out['perc2'][x] = m out['perc98'][x] = M else: out['perc2'][g] = m out['perc98'][g] = M out['X'][:,g] = (tmp - m)/(M - m) for cf in class_field: out[cf] = self.ref_db[cf].to_numpy(dtype=int) out['groups'] = self.ref_db['polygon_id'].to_numpy(dtype=int) return out def tiled_data(self, normalize=None): 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.int32) # Probably sub-optimal and error-prone, uses mean value from tile only imp = SimpleImputer(missing_values=np.nan, strategy='mean') X = imp.fit_transform(tile_data[vars].to_numpy()) if normalize is not None: X = (X - normalize[0]) / (normalize[1] - normalize[0]) yield tilenum,L,X def create_new_map(self): self.output_maps.append(np.zeros((self.H, self.W), dtype=np.int)) return 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) 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 if output_file is None: self.output_map[-1][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 not os.path.exists(os.path.dirname(output_file)): os.makedirs(os.path.dirname(output_file)) 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.uint16) as dst: # issue with compress option? #transform=self.transform, count=1, dtype=np.uint16, 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): 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 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, class_field='class'): 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(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] return true[pred>0], pred[pred>0] def tiled_objects(self, on_ref=False): assert(self.tiles is not None) idx = -1 if on_ref else 0 r = otb.itkRegion() for tn, t in self.tiles.items(): r['index'][0], r['index'][1] = t[0], t[1] r['size'][0], r['size'][1] = t[2], t[3] self.ref_obj_layer_pipe[idx].PropagateRequestedRegion('out', r) arr = self.ref_obj_layer_pipe[idx].GetImageAsNumpyArray('out').astype(np.uint32) yield tn, arr