"""
Things to consider adding:
choice of las or laz output
set default values for lasground_new params
clip structures step
lasclassify params to identify buildings
use veg polygon (if given) instead of inverse ground polygon to clip veg points
"""
from .lastools_fun import *
[docs]class DF(pd.DataFrame):
"""Extended pandas DataFrame class with an additional title attribute"""
def __init__(self, data=None, index=None, columns=None, dtype=None, copy=False, title=None):
pd.DataFrame.__init__(self, data, index, columns, dtype, copy)
self.title = title
def show(self):
if self.title:
print(self.title)
[docs]def ar1_acorr(series, maxlags=''):
"""Returns lag, autocorrelation, and confidence interval using geometric autocorrelation for AR1 fit of series"""
n = len(series)
if maxlags == '':
maxlags = int(n/2)
# use phi as lag-1 correlation of data
phi = np.corrcoef(series[:-1], series[1:])[0][1]
lags = range(maxlags+1)
acorrs = [phi**k for k in lags]
lower_band, upper_band = zip(*[r_confidence_interval(phi**k, n - k) for k in lags])
return lags, acorrs, lower_band, upper_band
[docs]def cmd(command):
"""Executes command prompt command"""
try:
res = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
except:
msg = 'Command failed: %s' % command
logging.error(msg)
raise Exception(msg)
msg = str(res.communicate()[1])
# if using for LAStools, get rid of the annoying LAStools licensing message.
msg = msg.replace(msg_from_lastools, '')
logging.info(msg)
return
[docs]def cox_acorr(series, maxlags=''):
"""
:param series: (list)
:param maxlags: (str)
:return: two lists (lags and autocorrelation), using Cox variant 3 of ACF
"""
n = len(series)
if maxlags == '':
maxlags = int(n/2)
xbar = np.mean(series)
lags = range(maxlags+1)
acorrs = []
for k in lags:
if k == 0:
acorrs.append(1)
else:
s1 = series[:-k]
s2 = series[k:]
numerator = 1.0/(n - k) * sum([(x1 - xbar) * (x2 - xbar) for x1, x2 in zip(s1, s2)])
denominator = 1.0/n * sum([(xi - xbar)**2 for xi in series])
acorrs.append(numerator*1.0/denominator)
return lags, acorrs
[docs]def ft(x, y):
"""Returns the fourier transform magnitude of the x,y data"""
n = len(x)
spacing = abs(x[1]-x[0])
xf = np.linspace(0, 1/(2.0*spacing), n//2)
yf = np.fft.fft(y)
yf = list(map(lambda k: 2.0/n*np.abs(k), yf))[:n//2]
return xf, yf
# input working directory for LAStools and directory containing .las/.laz files
# creates a .txt file for LAStools containing list of .las/.laz file names
# returns the name of the .txt file.
[docs]def lof_text(pwd, src):
"""creates a .txt file in pwd (LAStools bin) containing a list of .las/.laz filenames from src directory"""
filename = pwd + 'file_list.txt'
f = open(filename, 'w+')
if type(src) == str:
for i in get_las_files(src):
f.write('%s\n' % i)
else:
# this is the case when there are multiple source folders
for i in [name for source in src for name in get_las_files(source)]:
f.write('%s\n' % i)
f.close()
return filename
# input .las/.laz filename, outputs point density (after running lasinfo)
[docs]def pd(filename):
"""returns point density from lasinfo output .txt file"""
# name of txt output file from lasinfo
filename = filename[:-4] + '.txt'
f = open(filename, 'r')
text = f.readlines()
for line in text:
if line.startswith('point density:'):
d = line.split(' ')
d = d[d.index('returns') + 1]
return float(d)
# the main function that runs when 'run' button is clicked
@err_info
def process_lidar(lastoolsdir,
lidardir,
ground_poly,
cores,
units_code,
keep_orig_pts,
coarse_step,
coarse_bulge,
coarse_spike,
coarse_down_spike,
coarse_offset,
fine_step,
fine_bulge,
fine_spike,
fine_down_spike,
fine_offset
):
"""Executes main LAStools processing workflow. See readme for more info."""
classes = ['01-Default',
'02-Ground',
'05-Vegetation',
'06-Building'
]
if (ground_poly != '') and (keep_orig_pts == True):
# run on coarse and fine settings, need to clip and remove duplicates after merging
outdirs = ['00_separated',
'00_declassified',
'01_tiled',
'02a_lasground_new_coarse',
'02b_lasground_new_fine',
'03a_lasheight_coarse',
'03b_lasheight_fine',
'04a_lasclassify_coarse',
'04b_lasclassify_fine',
'05a_lastile_rm_buffer_coarse',
'05b_lastile_rm_buffer_fine',
'06a_separated_coarse',
'06b_separated_fine',
'07a_ground_clipped_coarse',
'07b_ground_clipped_fine',
'08_ground_merged',
'09_ground_rm_duplicates',
'10_veg_new_merged',
'11_veg_new_clipped',
'12_veg_merged',
'13_veg_rm_duplicates'
]
elif (ground_poly == '') and keep_orig_pts:
# only classify with coarse settings, no clipping, but need to remove duplicates
outdirs = ['00_separated',
'00_declassified',
'01_tiled',
'02_lasground_new',
'03_lasheight',
'04_lasclassify',
'05_lastile_rm_buffer',
'06_separated',
'08_ground_merged',
'09_ground_rm_duplicates',
'12_veg_merged',
'13_veg_rm_duplicates'
]
elif (ground_poly == '') and not keep_orig_pts:
# only classify with coarse setting, no clipping or removing duplicates necessary
outdirs = ['00_separated',
'00_declassified',
'01_tiled',
'02_lasground_new',
'03_lasheight',
'04_lasclassify',
'05_lastile_rm_buffer',
'06_separated'
]
elif (ground_poly != '') and not keep_orig_pts:
# run on coarse and fine settings, clip, but no removing duplicates needed
outdirs = ['00_separated',
'00_declassified',
'01_tiled',
'02a_lasground_new_coarse',
'02b_lasground_new_fine',
'03a_lasheight_coarse',
'03b_lasheight_fine',
'04a_lasclassify_coarse',
'04b_lasclassify_fine',
'05a_lastile_rm_buffer_coarse',
'05b_lastile_rm_buffer_fine',
'06a_separated_coarse',
'06b_separated_fine',
'07a_ground_clipped_coarse',
'07b_ground_clipped_fine',
'08_ground_merged',
'10_veg_new_merged',
'11_veg_new_clipped',
'12_veg_merged'
]
# make new directories for output from each step in processing
for outdir in outdirs:
if os.path.isdir(lidardir + outdir) == False:
os.mkdir(lidardir + outdir)
if len(os.listdir(lidardir + outdirs[0])) != 0:
msg = 'Output directories must initially be empty. Move or delete the data currently in output directories.'
logging.error(msg)
raise Exception(msg)
# in each 'separated' folder, create subdirs for each class type
if ground_poly != '':
sepdirs = [lidardir + '00_separated',
lidardir + '06a_separated_coarse',
lidardir + '06b_separated_fine'
]
else:
sepdirs = [lidardir + '00_separated',
lidardir + '06_separated'
]
for sepdir in sepdirs:
for class_type in classes:
class_dir = sepdir + '/' + class_type
if os.path.isdir(class_dir) == False:
os.mkdir(class_dir)
logging.info('Created directories for output data')
# create declassified points
logging.info('Declassifying copy of original point cloud...')
# get list of filenames for original LiDAR data (all .las and .laz files in lidardir)
lidar_files = []
for path, subdirs, files in os.walk(lidardir):
for name in files:
if name.endswith('.las') or name.endswith('.laz'):
lidar_files.append(path + '/' + name)
if lidar_files == []:
msg = 'No .las or .laz files in %s or its subdirectories' % lidardir
logging.error(msg)
raise Exception(msg)
# copy original files into '00_declassified' folder
for name in lidar_files:
shutil.copyfile(name, lidardir + '00_declassified/' + os.path.basename(name))
# make list of files for LASTools to process
lof = lof_text(lastoolsdir, lidardir + '00_declassified/')
# call LAStools command to declassify points and get point density
cmd('%slasinfo.exe -lof %s -set_classification 1 -otxt -cd' % (lastoolsdir, lof))
logging.info('OK')
# separate original data by class type
logging.info('Separating original data by class type...')
filename = lastoolsdir + 'file_list.txt'
f = open(filename, 'w+')
for i in lidar_files:
f.write('%s\n' % i)
f.close()
lof = filename
for class_type in classes:
odir = lidardir + '00_separated' + '/' + class_type + '/'
class_code = int(class_type.split('-')[0])
cmd('%slas2las.exe -lof %s -cores %i -keep_classification %i -odir %s -olas' % (
lastoolsdir, lof, cores, class_code, odir))
logging.info('OK')
# create tiling (max 1.5M pts per tile)
logging.info('Creating tiling...')
# get point density for each .las file
ds = []
for filename in get_las_files(lidardir + '00_declassified/'):
ds.append(pd(filename))
# use max point density out of all files to determine tile size
max_d = max(ds)
# width of square tile so we have max of 1.5M pts per tile (assuming same number of points per tile)
# throw in another factor of 0.5 to make sure tiles will be small enough, round to nearest 10
tile_size = round(0.5 * np.sqrt((1.5 * 10 ** 6) / max_d), -1)
logging.info('Using tile size of %i' % tile_size)
odir = lidardir + '01_tiled/'
# call LAStools command to create tiling
cmd('%slastile.exe -lof %s -cores %i -o tile.las -tile_size %i -buffer 5 -faf -odir %s -olas' % (
lastoolsdir, lof, cores, tile_size, odir))
logging.info('OK')
# check to make sure tiles are small enough
logging.info('Checking if largest file has < 1.5M pts (to avoid licensing restrictions)...')
largest_file = get_largest(odir)
num = pts(largest_file, lastoolsdir)
if num < 1500000:
logging.info('Largest file has %i points, tiles small enough.' % num)
else:
logging.info('Tile size not small enough. Retrying with a smaller tile size...')
while num >= 1500000:
# delete original set of tiles
folder = odir
for the_file in os.listdir(folder):
file_path = os.path.join(folder, the_file)
try:
if os.path.isfile(file_path):
os.unlink(file_path)
except:
logging.warning('Couldn\'nt delete %s' % file_path)
# redo tiling
tile_size = int(tile_size * num * 1.0 / 1500000)
logging.info('Using tile size of %i' % tile_size)
cmd('%slastile.exe -lof %s -cores %i -o tile.las -tile_size %i -buffer 5 -faf -odir %s -olas' % (
lastoolsdir, lof, cores, tile_size, odir))
# recheck largest tile number of points
logging.info('Checking if largest file has < 1.5M pts (to avoid licensing restrictions)...')
largest_file = get_largest(odir)
num = pts(largest_file, lastoolsdir)
if num >= 1500000:
logging.info('Tile size not small enough. Retrying with a smaller tile size...')
logging.info('OK')
# run lasground_new on coarse and fine settings
logging.info('Running ground classification on coarse setting...')
lof = lof_text(lastoolsdir, lidardir + '01_tiled/')
if ground_poly != '':
odir = lidardir + '02a_lasground_new_coarse/'
else:
odir = lidardir + '02_lasground_new'
cmd(
'%slasground_new.exe -lof %s -cores %i %s -step %s -bulge %s -spike %s -down_spike %s -offset %s -hyper_fine -odir %s -olas' % (
lastoolsdir,
lof,
cores,
units_code,
coarse_step,
coarse_bulge,
coarse_spike,
coarse_down_spike,
coarse_offset,
odir
)
)
logging.info('OK')
if ground_poly != '':
logging.info('Running ground classification on fine setting...')
odir = lidardir + '02b_lasground_new_fine/'
cmd(
'%slasground_new.exe -lof %s -cores %i %s -step %s -bulge %s -spike %s -down_spike %s -offset %s -hyper_fine -odir %s -olas' % (
lastoolsdir,
lof,
cores,
units_code,
fine_step,
fine_bulge,
fine_spike,
fine_down_spike,
fine_offset,
odir
)
)
logging.info('OK')
# run lasheight on each data set
logging.info('Measuring height above ground for non-ground points...')
if ground_poly != '':
lof = lof_text(lastoolsdir, lidardir + '02a_lasground_new_coarse/')
odir = lidardir + '03a_lasheight_coarse/'
else:
lof = lof_text(lastoolsdir, lidardir + '02_lasground_new/')
odir = lidardir + '03_lasheight/'
cmd('%slasheight.exe -lof %s -cores %i -odir %s -olas' % (lastoolsdir, lof, cores, odir))
if ground_poly != '':
lof = lof_text(lastoolsdir, lidardir + '02b_lasground_new_fine/')
odir = lidardir + '03b_lasheight_fine/'
cmd('%slasheight.exe -lof %s -cores %i -odir %s -olas' % (lastoolsdir, lof, cores, odir))
logging.info('OK')
# run lasclassify on each data set
logging.info('Classifying non-ground points on coarse setting...')
if ground_poly != '':
lof = lof_text(lastoolsdir, lidardir + '03a_lasheight_coarse/')
odir = lidardir + '04a_lasclassify_coarse/'
else:
lof = lof_text(lastoolsdir, lidardir + '03_lasheight/')
odir = lidardir + '04_lasclassify/'
cmd('%slasclassify.exe -lof %s -cores %i %s -odir %s -olas' % (lastoolsdir, lof, cores, units_code, odir))
logging.info('OK')
if ground_poly != '':
logging.info('Classifying non-ground points on fine setting...')
lof = lof_text(lastoolsdir, lidardir + '03b_lasheight_fine/')
odir = lidardir + '04b_lasclassify_fine/'
cmd('%slasclassify.exe -lof %s -cores %i %s -odir %s -olas' % (lastoolsdir, lof, cores, units_code, odir))
logging.info('OK')
# remove tile buffers on each data set
logging.info('Removing tile buffers...')
if ground_poly != '':
lof = lof_text(lastoolsdir, lidardir + '04a_lasclassify_coarse/')
odir = lidardir + '05a_lastile_rm_buffer_coarse/'
else:
lof = lof_text(lastoolsdir, lidardir + '04_lasclassify/')
odir = lidardir + '05_lastile_rm_buffer/'
cmd('%slastile.exe -lof %s -cores %i -remove_buffer -odir %s -olas' % (lastoolsdir, lof, cores, odir))
if ground_poly != '':
lof = lof_text(lastoolsdir, lidardir + '04b_lasclassify_fine/')
odir = lidardir + '05b_lastile_rm_buffer_fine/'
cmd('%slastile.exe -lof %s -cores %i -remove_buffer -odir %s -olas' % (lastoolsdir, lof, cores, odir))
logging.info('OK')
# separate into files for each class type
logging.info('Separating points by class type on coarse setting...')
# coarse
if ground_poly != '':
lof = lof_text(lastoolsdir, lidardir + '05a_lastile_rm_buffer_coarse/')
podir = lidardir + '06a_separated_coarse'
else:
lof = lof_text(lastoolsdir, lidardir + '05_lastile_rm_buffer/')
podir = lidardir + '06_separated'
for class_type in classes:
odir = podir + '/' + class_type + '/'
class_code = int(class_type.split('-')[0])
cmd('%slas2las.exe -lof %s -cores %i -keep_classification %i -odir %s -olas' % (
lastoolsdir, lof, cores, class_code, odir))
logging.info('OK')
if ground_poly == '' and (keep_orig_pts == False):
ground_results = podir + '/' + '02-Ground' + '/'
veg_results = podir + '/' + '05-Vegetation' + '/'
if ground_poly != '':
logging.info('Separating points by class type on fine setting...')
# fine
lof = lof_text(lastoolsdir, lidardir + '05b_lastile_rm_buffer_fine/')
for class_type in classes:
odir = lidardir + '06b_separated_fine' + '/' + class_type + '/'
class_code = int(class_type.split('-')[0])
cmd('%slas2las.exe -lof %s -cores %i -keep_classification %i -odir %s -olas' % (
lastoolsdir, lof, cores, class_code, odir))
logging.info('OK')
# clip ground data sets with ground polygon
if ground_poly != '':
logging.info('Clipping ground points to inverse ground polygon on coarse setting...')
# keep points outside ground polygon for coarse setting (-interior flag)
lof = lof_text(lastoolsdir, lidardir + '06a_separated_coarse' + '/' + '02-Ground' + '/')
odir = lidardir + '07a_ground_clipped_coarse/'
cmd('%slasclip.exe -lof %s -cores %i -poly %s -interior -donuts -odir %s -olas' % (
lastoolsdir, lof, cores, ground_poly, odir))
logging.info('OK')
logging.info('Clipping ground points to ground polygon on fine setting...')
# keep points inside ground polygon for fine setting
lof = lof_text(lastoolsdir, lidardir + '06b_separated_fine' + '/' + '02-Ground' + '/')
odir = lidardir + '07b_ground_clipped_fine/'
cmd('%slasclip.exe -lof %s -cores %i -poly %s -donuts -odir %s -olas' % (
lastoolsdir, lof, cores, ground_poly, odir))
logging.info('OK')
# merge processed ground points with original data set ground points
if keep_orig_pts:
logging.info('Merging new and original ground points...')
if ground_poly != '':
sources = [lidardir + '07a_ground_clipped_coarse/', lidardir + '07b_ground_clipped_fine/',
lidardir + '00_separated' + '/' + '02-Ground' + '/']
else:
sources = [lidardir + '06_separated' + '/' + '02-Ground' + '/',
lidardir + '00_separated' + '/' + '02-Ground' + '/']
# just use new points
elif ground_poly != '':
logging.info('Merging new ground points...')
sources = [lidardir + '07a_ground_clipped_coarse/', lidardir + '07b_ground_clipped_fine/']
if keep_orig_pts or (ground_poly != ''):
lof = lof_text(lastoolsdir, sources)
odir = lidardir + '08_ground_merged/'
ground_results = odir # will be overwritten if rm_duplicates block runs
cmd('%slastile.exe -lof %s -cores %i -o tile.las -tile_size %i -faf -odir %s -olas' % (
lastoolsdir, lof, cores, tile_size, odir))
logging.info('OK')
# remove duplicate ground points
if keep_orig_pts:
logging.info('Removing duplicate ground points...')
lof = lof_text(lastoolsdir, lidardir + '08_ground_merged/')
odir = lidardir + '09_ground_rm_duplicates/'
ground_results = odir
cmd('%slasduplicate.exe -lof %s -cores %i -lowest_z -odir %s -olas' % (lastoolsdir, lof, cores, odir))
logging.info('OK')
# merge new veg points
if ground_poly != '':
logging.info('Merging new vegetation points from coarse and fine run...')
sources = [lidardir + '06a_separated_coarse' + '/' + '05-Vegetation' + '/',
lidardir + '06b_separated_fine' + '/' + '05-Vegetation' + '/']
lof = lof_text(lastoolsdir, sources)
odir = lidardir + '10_veg_new_merged/'
cmd('%slastile.exe -lof %s -cores %i -o tile.las -tile_size %i -faf -odir %s -olas' % (
lastoolsdir, lof, cores, tile_size, odir))
logging.info('OK')
# clip new veg points
# keeping points outside the ground polygon
logging.info('Clipping new vegetation points...')
lof = lof_text(lastoolsdir, lidardir + '10_veg_new_merged/')
odir = lidardir + '11_veg_new_clipped/'
cmd('%slasclip.exe -lof %s -cores %i -poly %s -interior -donuts -odir %s -olas' % (
lastoolsdir, lof, cores, ground_poly, odir))
logging.info('OK')
# merge with original veg points
if keep_orig_pts:
logging.info('Merging new and original vegetation points...')
if ground_poly != '':
sources = [lidardir + '11_veg_new_clipped/', lidardir + '00_separated' + '/' + '05-Vegetation' + '/']
else:
sources = [lidardir + '06_separated' + '/' + '05-Vegetation' + '/',
lidardir + '00_separated' + '/' + '05-Vegetation' + '/']
elif ground_poly != '':
logging.info('Retiling new vegetation points...')
sources = [lidardir + '11_veg_new_clipped/']
if (keep_orig_pts == True) or (ground_poly != ''):
lof = lof_text(lastoolsdir, sources)
odir = lidardir + '12_veg_merged/'
veg_results = odir # will be overwritten if rm_duplicates block runs
cmd('%slastile.exe -lof %s -cores %i -o tile.las -tile_size %i -faf -odir %s -olas' % (
lastoolsdir, lof, cores, tile_size, odir))
logging.info('OK')
# remove duplicate veg points
if keep_orig_pts:
logging.info('Removing duplicate vegetation points...')
lof = lof_text(lastoolsdir, lidardir + '12_veg_merged/')
odir = lidardir + '13_veg_rm_duplicates/'
veg_results = odir
cmd('%slasduplicate.exe -lof %s -cores %i -lowest_z -odir %s -olas' % (lastoolsdir, lof, cores, odir))
logging.info('OK')
logging.info('Processing finished.')
logging.info('Outputs in:')
logging.info('%s\n%s' % (ground_results, veg_results))
return
[docs]def pts(filename, lastoolsdir):
"""returns number of points in las file"""
# call lasinfo on the file
cmd('%slasinfo.exe -i %s -otxt -histo number_of_returns 1' % (lastoolsdir, filename))
# name of txt output file from lasinfo
txt = filename[:-4] + '.txt'
f = open(txt, 'r')
text = f.readlines()
for line in text:
if line.endswith('element(s)\n'):
d = line.split(' ')
d = d[d.index('for') + 1]
return int(d)
def r_to_z(r):
return 0.5 * np.log((1+r)*1.0/(1-r))
[docs]def r_confidence_interval(r, n, alpha=0.05):
"""Retrieves the confidence interval at the 1-alpha level for correlation of r with n observations
when alpha=0.05, it returns the range of possible population correlations at the 95% confidence level
so if 0 is not within the bounds, then the correlation is statistically significant at the 95% level
:param r: correlation (float)
:param n: number of observations (int)
:param alpha: confidence level (float)
:return: Confidence interval (low and high) as sequence (list or tuple) of floats.
"""
if r == 1:
return 1, 1
z = r_to_z(r)
se = 1.0 / np.sqrt(n - 3)
z_crit = stats.norm.ppf(1 - alpha / 2) # 2-tailed z critical value
lo = z - z_crit * se
hi = z + z_crit * se
return z_to_r(lo), z_to_r(hi)
def white_noise_confidence_interval(n):
return -1.0/n - 2.0/np.sqrt(n), -1.0/n + 2.0/np.sqrt(n)
[docs]def white_noise_acf_ci(series, maxlags=''):
"""Returns the 95% confidence interval for white noise ACF"""
n = len(series)
if maxlags == '':
maxlags = int(n/2)
lags = range(maxlags+1)
lims = [white_noise_confidence_interval(n-k) for k in lags]
lower_lims, upper_lims = list(zip(*lims))
return lags, lower_lims, upper_lims
def z_to_r(z):
return (np.exp(2*z)-1)*1.0/(np.exp(2*z) + 1)