Examples

Process module

The submodule process allows full control over all types of workflows through use of a number of processing flags. A minimum working example usage of the the PySESA module, accepting all default values for parameters, is:

import pysesa
infile = `/home/me/mypointcloudfile.txt'
pysesa.process(infile)

This instance writes out the following results file whose name contains some of the processing parameters:

/home/me/mypointcloudfile.txt_zstat_detrend4_outres0.5_proctype1_mxpts512_minpts16.xyz

Test module

The above is the same as passing a list of default-valued variables to process, which is included for completeness in the test module:

out = 1               # 1 m output grid
detrend = 4           # detrend type: ODR plane
proctype = 1          # Processing type: spectral parameters (no smoothing) only
mxpts = 1024          # Maximum points per window
res = 0.05            # 5 cm grid resolution for detrending and spectral analysis
nbin = 20             # Number of bins for spectral binning
lentype = 1           # Integral lengthscale type: l<0.5
taper = 1             # Hann taper before spectral analysis
prc_overlap = 0       # No overlap between successive windows
minpts = 64           # Minimum points per window

pysesa.process(infile, out, detrend, proctype, mxpts, res, nbin, lentype, minpts, taper, prc_overlap)

Minimal example on 1 window

A minimal example analysis of spatial and spectral analysis on just 1 window of data:

# import module
import pysesa

# read point cloud from file
pointcloud = pysesa.read.txtread(infile)

# create windows of data
windows = pysesa.partition(pointcloud).getdata()

# process window number 50
k=50

# get all spectral statistics for that window
spec_stats = pysesa.spectral(pointcloud[windows[k],:3].astype('float64')).getdata()

# get all spatial statistics for that window
spat_stats = pysesa.spatial(pointcloud[windows[k],:3].astype('float64')).getdata()

Minimal example, all windows using parallel processing

and to extend this to all windows, utilising parallel processing over all available cores, could be achieved using the following minimal example:

# define a function that will get repeatedly read by the parallel processing queue
def get_spat_n_spec(pts):
   return pysesa.spatial(pts.astype('float64')).getdata() + pysesa.spectral(pts.astype('float64')).getdata()

# import the parallel processing libraries
from joblib import Parallel, delayed, cpu_count

# Processing type: spatial plus spectral parameters (no smoothing)
proctype = 4

# process each window with all available cores, by queueing each window in a sequence
# and processing until they are all done
w = Parallel(n_jobs=cpu_count(), verbose=0)(delayed(get_spat_n_spec)(pointcloud[windows[k],:3])
      for k in xrange(len(windows)))

# parse out the outputs into variables
x, y, z_mean, z_max, z_min, z_range, sigma, skewness, kurtosis, n, ...
slope, intercept, r_value, p_value, std_err, d, l, wmax, wmean, rms1, rms2, ...
Z, E, sigma, T0_1, T0_2, sw1, sw2, m0, m1, m2, m3, m4, phi = zip(*w)

Lengthscale module

To obtain just the integral lengthscale of the kth window, detrended using the orthogonal distance regression detrending technique, one could use:

detrend = 4   # Orthogonal distance regression
pysesa.lengthscale(pysesa.detrend(pointcloud[windows[k],:3],detrend).getdata()).getlengthscale()

Spatial module

and to get the spatial statistics from the same data:

pysesa.spatial(pysesa.detrend(pointcloud[windows[k],:3],detrend).getdata()).getdata()

Spectral module

Here, the output grid resolution is changed to 25 cm, and the various outputs from the spectral module are obtained separately:

# 25 cm output grid
out = 0.25

# re-create windows of data
windows = pysesa.partition(pointcloud, out).getdata()

result = pysesa.spectral(pointcloud[windows[k],:3].astype('float64'))

# get all spectral parameters
result.getdata()

# get the fit parameters for log-log power spectrum
result.getpsdparams()

# get integral lengthscale
result.getlengthscale()

# get spectral moment parameters
result.getmoments()

# get rms and wavelength parameters
result.getlengths()

Plot module

This assumes you have run the process module and have an output file (‘/home/my_pysesa_output_file.xyz’):

# load pysesa
import pysesa

# create a pysesa::plot instance
p = pysesa.plot('/home/my_pysesa_output_file.xyz')

# create a 3d plot of the point cloud
p.plt_xyz()

# create a 2d plot of the gridded surface from the decimated point cloud
p.grd_xyz()

# create a 3d plot of the gridded surface from the decimated point cloud
# colour-coded by amplitude
p.grd_xyz3d()

# create a 3d plot of the decimated spectral slope
# colour-coded by amplitude
p.plt_xy_var('slope')

# create a 3d plot of all output decimated parameters
# colour-coded by amplitude
p.plt_xy_vars()

# create a 3d plot of decimated fractal dimension
# gridded and colour-coded by amplitude
p.grd_var_3d('d')

# create a 3d plot of all output decimated parameters
# gridded and colour-coded by amplitude
p.grd_vars_3d()

# plot also supports data retrieval
# retrieve the original point cloud
xyz = p.get_xyz()

# retrieve the decimated point cloud of all parameters
pc = p.get_pc()

# retrieve the decimated point cloud of all parameters
# in dictionary format
pc_dict = p.parse_pc_vars()
# then show what's in there
print pc.dict.keys()