florent/ld2daps
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Better filtering and cliping when loading lidar

Browse Source
This commit is contained in:
Florent Guiotte 2018-08-29 17:20:54 +02:00
parent eb6a39896c
commit 11baca069b
2 changed files with 25 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Notebooks/Raster Factory.ipynb
View File

@ -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')"
]
},
{

43
raster_assistant.py
View File

@ -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):