diff --git a/Notebooks/Raster Factory.ipynb b/Notebooks/Raster Factory.ipynb index 4596539..25ddbc4 100644 --- a/Notebooks/Raster Factory.ipynb +++ b/Notebooks/Raster Factory.ipynb @@ -40,7 +40,7 @@ "metadata": {}, "outputs": [], "source": [ - "C1 = ra.bulk_load('../Data/lidar/C1', 'C1', filter_treshold=0, dtype=np.float32)\n", + "C1 = ra.bulk_load('../Data/lidar/C1', 'C1', filter_treshold=.01, clip_treshold=.1, dtype=np.float32)\n", "#C2 = ra.bulk_load('../Data/lidar/C2', 'C2', filter_treshold=.5, dtype=np.float32)\n", "#C3 = ra.bulk_load('../Data/lidar/C3', 'C3', filter_treshold=.5, dtype=np.float32)\n", "#C123 = ra.bulk_load('../Data/lidar', 'C123', filter_treshold=.5, dtype=np.float32)" @@ -59,7 +59,7 @@ "metadata": {}, "outputs": [], "source": [ - "raster = ra.rasterize_cache(C1, 'intensity', 1., 'nearest', True, cache_dir='../Res/')" + "raster = ra.rasterize_cache('intensity', C1, 1., 'nearest', False, cache_dir='../Res/enrichment_rasters')" ] }, { diff --git a/raster_assistant.py b/raster_assistant.py index c5d5a3b..25594ec 100644 --- a/raster_assistant.py +++ b/raster_assistant.py @@ -17,12 +17,12 @@ sys.path.append('../triskele/python/') import triskele import rasterizer -def rasterize_cache(data, field, resolution=1., method='linear', reverse_alt=False, cache_dir='/tmp'): +def rasterize_cache(field, data, resolution=.5, interpolation='nearest', reverse_direction=False, cache_dir='/tmp'): """Cache layer for rasterize""" cache_dir = Path(cache_dir) - name = '{}_{}_{}_{}{}'.format(data.name, field, str(resolution).replace('.', '_'), method, - '_reversed' if reverse_alt else '') + name = '{}_{}_f{}_c{}_r{}_{}{}'.format(field, data.name, data.filter_treshold, data.clip_treshold, resolution, + interpolation, '_reversed' if reverse_direction else '') png_file = cache_dir.joinpath(Path(name + '.png')) tif_file = cache_dir.joinpath(Path(name + '.tif')) @@ -30,7 +30,7 @@ def rasterize_cache(data, field, resolution=1., method='linear', reverse_alt=Fal print ('WARNING: Loading cached result {}'.format(tif_file)) raster = triskele.read(tif_file) else: - raster = rasterizer.rasterize(data.spatial, getattr(data, field), resolution, method, reverse_alt, np.float32) + raster = rasterizer.rasterize(data.spatial, getattr(data, field), resolution, interpolation, reverse_direction, np.float32) triskele.write(tif_file, raster) plt.imsave(png_file, raster) @@ -50,17 +50,20 @@ def find_las(path): return las def extremum_filter(data, treshold=.5): + """Return boolean filter of extremums in data, treshold in percentil [0;100]""" tresholds = np.percentile(data, [treshold, 100 - treshold]) return np.logical_or(data < tresholds[0], data > tresholds[1]) -def auto_filter(data, treshold=.5): +def extremum_clip(data, treshold=.5): + """Clip extremums in data, treshold in percentil [0;100]""" tresholds = np.percentile(data, [treshold, 100 - treshold]) - - data[data < tresholds[0]] = tresholds[0] - data[data > tresholds[1]] = tresholds[1] + np.clip(data, *tresholds, out=data) -def bulk_load(path, name=None, filter_treshold=.1, dtype=None): - data = {'file': path, 'name': name} +def bulk_load(path, name=None, filter_treshold=.1, clip_treshold=.5, dtype=None): + data = {'file': path, + 'name': name, + 'filter_treshold': filter_treshold, + 'clip_treshold': clip_treshold} attributes = ['x', 'y', 'z', 'intensity', 'num_returns']#, 'scan_angle_rank', 'pt_src_id', 'gps_time'] for a in attributes: @@ -84,19 +87,19 @@ def bulk_load(path, name=None, filter_treshold=.1, dtype=None): print('\rCreate matrices: [Done]') print('Filter data...', end='') - Z = data['z'].copy() - auto_filter(Z) - data['spatial'] = np.array((data['x'], data['y'], Z)).T - del Z - - t = .1 - efilter = np.logical_or(extremum_filter(data['intensity'], t), extremum_filter(data['z'])) - - attributes.append('spatial') + z_or_i = np.logical_or(extremum_filter(data['z'], filter_treshold), + extremum_filter(data['intensity'], filter_treshold)) for i, a in enumerate(attributes): print('\rFilter data: [{:3d}%]'.format(int(i/len(attributes) * 100)), end='') - data[a] = data[a][np.logical_not(efilter)] + data[a] = data[a][~z_or_i] print('\rFilter data: [Done]') + + print('Clip data...', end='') + extremum_clip(data['z'], clip_treshold) + extremum_clip(data['intensity'], clip_treshold) + data['spatial'] = np.array((data['x'], data['y'], data['z'])).T + attributes.append('spatial') + print('\rClip data: [Done]') class TMPLAS(object): def __init__(self, d):