105 lines
3.3 KiB
Python
105 lines
3.3 KiB
Python
#!/usr/bin/python
|
|
# -*- coding: utf-8 -*-
|
|
# \file raster_assistant.py
|
|
# \brief TODO
|
|
# \author Florent Guiotte <florent.guiotte@gmail.com>
|
|
# \version 0.1
|
|
# \date 03 juil. 2018
|
|
#
|
|
# TODO details
|
|
|
|
import sys
|
|
from pathlib import Path
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
import laspy
|
|
sys.path.append('../triskele/python/')
|
|
import triskele
|
|
import rasterizer
|
|
|
|
def rasterize_cache(data, field, resolution=1., method='linear', reverse_alt=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 '')
|
|
png_file = cache_dir.joinpath(Path(name + '.png'))
|
|
tif_file = cache_dir.joinpath(Path(name + '.tif'))
|
|
|
|
if tif_file.exists() :
|
|
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)
|
|
triskele.write(tif_file, raster)
|
|
plt.imsave(png_file, raster)
|
|
|
|
return raster
|
|
|
|
def find_las(path):
|
|
path = Path(path)
|
|
las = list()
|
|
|
|
if path.is_dir():
|
|
for child in path.iterdir():
|
|
las.extend(find_las(child))
|
|
|
|
if path.is_file() and path.suffix == '.las':
|
|
las.append(path)
|
|
|
|
return las
|
|
|
|
def extremum_filter(data, treshold=.5):
|
|
tresholds = np.percentile(data, [treshold, 100 - treshold])
|
|
return np.logical_or(data < tresholds[0], data > tresholds[1])
|
|
|
|
def auto_filter(data, treshold=.5):
|
|
tresholds = np.percentile(data, [treshold, 100 - treshold])
|
|
|
|
data[data < tresholds[0]] = tresholds[0]
|
|
data[data > tresholds[1]] = tresholds[1]
|
|
|
|
def bulk_load(path, name=None, filter_treshold=.1, dtype=None):
|
|
data = {'file': path, 'name': name}
|
|
|
|
attributes = ['x', 'y', 'z', 'intensity', 'num_returns']#, 'scan_angle_rank', 'pt_src_id', 'gps_time']
|
|
for a in attributes:
|
|
data[a] = list()
|
|
|
|
print('Load data...')
|
|
for f in find_las(path):
|
|
print('{}: '.format(f), end='')
|
|
infile = laspy.file.File(f)
|
|
for i, a in enumerate(attributes):
|
|
print('\r {}: [{:3d}%]'.format(f, int(i/len(attributes) * 100)), end='')
|
|
data[a].extend(getattr(infile, a))
|
|
infile.close()
|
|
print('\r {}: [Done]'.format(f))
|
|
|
|
|
|
print('Create matrices...', end='')
|
|
for i, a in enumerate(attributes):
|
|
print('\rCreate matrices: [{:3d}%]'.format(int(i/len(attributes) * 100)), end='')
|
|
data[a] = np.array(data[a], dtype=dtype)
|
|
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')
|
|
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)]
|
|
print('\rFilter data: [Done]')
|
|
|
|
class TMPLAS(object):
|
|
def __init__(self, d):
|
|
self.__dict__.update(d)
|
|
|
|
return TMPLAS(data) |