In [1]:
# FULL SCRIPT - SOME STUFF HERE WAS NOT USED IN THE FINAL PAPER
# LIKE R3.VOLUME AND THE NORMALITY TESTS OF DISPLACEMENT LENGTH
# 
##############################################################################
# Script used in the paper:
# Dune migration and volume change from airborne LiDAR, terrestrial LiDAR 
# and Structure from Motion--Multi View Stereo
# by
# Carlos H. Grohmann et al - 2019/2020
# guano (at) usp (dot) br
# Institute of Energy and Environment - University of Sao Paulo
#
# Please check the GitHub repo for the final reference to the paper
##############################################################################

In [2]:
# external file 'azim.py' has the function to calculate azimuth and length of lines

In [3]:
# Shapefiles provided with this script:
# area_dem - mask for ALS and SfM interpolation and comparison
# area_volume - mask for volume calculation (ALS and SfM)
# area_dunefield_als - mask for the entire dunefield (ALS coverage)
# area_tls - mask of TLS survey area
# area_tls_sfm - mask used to constrain random points (TLS x SfM)
# crests_als_2010 - crestlines of 2010 (ALS)
# crests_sfm_2019 - crestlines of 2019 (SfM)
# displacement_crests_2010_2019 - lines connecting 2010 and 2019 crests

In [4]:
# Point cloud data used in this paper is available via OpenTopography
# Airborne LiDAR (2010) - http://dx.doi.org/10.5069/G9DN430Z
# Sfm-MVS (2019) - http://dx.doi.org/10.5069/G9DV1H19
# Terrestrial LiDAR (2019) -

In [5]:
# Note on LiDAR data used in this paper:
# my GRASS wasn't working witb LibLAS for direct LAS import, so I converted to TXT.
# depending on the GRASS installation, you can use v.in.lidar or v.in.pdal as well

In [6]:
# import python libraries
import sys, os, itertools
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import matplotlib as mpl
import pandas as pd
import seaborn as sns
import subprocess
from IPython.display import Image # can use this to display GRASS maps
# stats
from statsmodels.graphics.gofplots import qqplot
from scipy.stats import linregress

In [7]:
# import module to calculate azimuth of lines -- requires python GDAL/OGR libraries 
# in my mac, I needed to set this or jupyter wouldn't find the ogr module:
sys.path.append('/usr/local/opt/gdal2-python/lib/python2.7/site-packages')
import azim

In [8]:
# helper func: round to nearest 5
def round5(x):
    rounded = int(round(x/5.0)*5.0)
    return rounded

In [9]:
# error measurements
def err_mse(x1, x2, axis=0):
    """mean squared error"""
    x1 = np.asanyarray(x1)
    x2 = np.asanyarray(x2)
    return np.mean((x1-x2)**2, axis=axis)

def err_rmse(x1, x2, axis=0):
    """root mean squared error"""
    x1 = np.asanyarray(x1)
    x2 = np.asanyarray(x2)
    return np.sqrt(err_mse(x1, x2, axis=axis))

def err_mae(x1, x2, axis=0):
    """mean absolute error"""
    x1 = np.asanyarray(x1)
    x2 = np.asanyarray(x2)
    return np.mean(np.abs(x1-x2), axis=axis)

In [10]:
def fdr(df):
    ''' Freedman–Diaconis rule to calculate optimal bin size of histograms
    '''
    q25 = df.quantile(0.25)
    q75 = df.quantile(0.75)
    len_max = df.max()
    len_min = df.min()
    len_n = len(df)
    IQR = q75 - q25
    bin_width = 2. * IQR * (len_n)**(-1./3.)
    nbins = (len_max - len_min)/bin_width
    return int(nbins)

In [11]:
# aux func
def raster_as_array(raster):
    ''' return GRASS raster as numpy array - remove null values '''
    grass.run_command('g.region', raster=raster, flags='a')
    raster_array = garray.array()
    raster_array.read(raster, null=np.nan)
    raster_array_flat = raster_array.flatten(order='C')
    raster_1d = raster_array_flat[~np.isnan(raster_array_flat)]
    return raster_1d

In [12]:
# matplotlib figures appear inline in the notebook rather than in a new window.
%matplotlib inline

In [13]:
# create GRASS GIS runtime environment
# with this, you can run GRASS without startig a shell/gui session
gisbase = subprocess.check_output(["grass76", "--config", "path"]).strip()
os.environ['GISBASE'] = gisbase
sys.path.append(os.path.join(gisbase, "etc", "python"))

# GRASS GIS imports
import grass.script as grass
import grass.script.setup as gsetup
import grass.script.array as garray
import grass.script.vector as gvect

In [14]:
# set GRASS GIS session data
# I use two systems, so this makes things a bit easier
if sys.platform == "linux" or sys.platform == "linux2":
    rcfile = gsetup.init(gisbase, "/mnt/sda/grassdata/", "utm", "garopaba_22J")
elif sys.platform == "darwin":
    rcfile = gsetup.init(gisbase, "/Volumes/MacintoshHD2/grassdata/", "utm", "garopaba_22J")
# elif platform == "win32":
    # Windows...
    
# grass.message('Current GRASS GIS 7 environment:')
# print grass.gisenv()

In [15]:
# overwrite for GRASS modules
ow = True

In [16]:
# check region
# grass.parse_command('g.region', flags='pg')
# grass.run_command('g.region', flags='pg')

In [17]:
# Data dir 
# use this to set different paths for different systems
if sys.platform == "linux" or sys.platform == "linux2":
    dataDir = '/mnt/sda/Dropbox/USP/projetosPesquisa/LiDAR_terrestre_SfM/_areas_estudo/garopaba/DEMs_shapes/'
elif sys.platform == "darwin":
    dataDir = '/Volumes/MacintoshHD2/Dropbox/USP/projetosPesquisa/LiDAR_terrestre_SfM/_areas_estudo/garopaba/DEMs_shapes/'
#     dataDir = '_path_to_your_files_'

In [18]:
###########################################################################
###########################################################################
# 
# Import vector masks, import point clouds as raster to create maps 
# of point density/m2
# 
###########################################################################
###########################################################################

In [19]:
# ------------------------------------------------------------------------
# vector masks
# ------------------------------------------------------------------------

In [20]:
# ALS - entire dune field (polygon, manually digitized)
shp_als_mask = dataDir+'area_garopaba.shp'
mask_als = 'mask_als'

In [21]:
grass.run_command('v.in.ogr', input=shp_als_mask, output=mask_als, flags='o', overwrite=ow)

In [22]:
# TLS - dune field TLS area mask (polygon, manually digitized)
# import points in GRASS with r.in.xyz with 1m resolution, export as tiff, draw polygon in QGIS:
# g.region -pa res=01 w=732187 e=732674 n=6900302 s=6899889
# r.in.xyz --overwrite input=pcloud_tls_20mm.txt output=mdt_tls_2019 method=mean separator=comma
# r.out.gdal --overwrite input=mdt_tls_2019@garopaba_22J output=pcloud_1m.tif format=GTiff
shp_tls_mask = dataDir+'area_tls.shp'
mask_tls = 'mask_tls'

In [23]:
grass.run_command('v.in.ogr', input=shp_tls_mask, output=mask_tls, flags='o', overwrite=ow)

In [24]:
# SfM - dune field area mask (polygon, manually digitized)
shp_sfm_mask = dataDir+'area_dem.shp'
mask_sfm = 'mask_sfm' # mask for clipping (pretty much the SfM area)

In [25]:
grass.run_command('v.in.ogr', input=shp_sfm_mask, output=mask_sfm, flags='o', overwrite=ow)

In [26]:
# ------------------------------------------------------------------------
# create maps of point density / m2, adjust colors (hist equalize)
# ------------------------------------------------------------------------

In [27]:
# ALS
pcloud_als_txt = '_path_to_your_files_/pcloud_als.txt' 
als_density = 'als_density'
grass.run_command('g.region', vector=mask_sfm, res=1, flags='pa')

In [28]:
grass.run_command('r.in.xyz', input=pcloud_als_txt, output=als_density, method='n', separator='comma', overwrite=ow) 
grass.run_command('r.null', map=als_density, setnull=0)
grass.run_command('r.colors', map=als_density, color='viridis', flags='e')

In [29]:
# display with virtual monitor rendered to png
grass.run_command('d.mon', start='png', output='view.png', overwrite=ow)
grass.run_command('d.rast', map=als_density)
grass.run_command('d.vect', map=mask_sfm, fill_color='none')
grass.run_command('d.legend', raster=als_density, at='5,80,0,4', flags='sf')
Image('view.png')


Out[29]:

In [30]:
grass.run_command('d.mon', stop='png')


Out[30]:
0

In [31]:
# SfM - full resolution
pcloud_sfm_full = '_path_to_your_files_/pcloud_sfm_full.txt'
sfm_full_density = 'sfm_full_density'
grass.run_command('g.region', vector=mask_sfm, res=1, flags='pa')

In [32]:
grass.run_command('r.in.xyz', input=pcloud_sfm_full, output=sfm_full_density, method='n', separator='comma', overwrite=ow)
grass.run_command('r.null', map=sfm_full_density, setnull=0)
grass.run_command('r.colors', map=sfm_full_density, color='viridis', flags='e')

In [33]:
# display with virtual monitor rendered to png
grass.run_command('d.mon', start='png', output='view.png', overwrite=ow)
grass.run_command('d.rast', map=sfm_full_density)
grass.run_command('d.vect', map=mask_sfm, fill_color='none')
grass.run_command('d.legend', raster=sfm_full_density, at='5,80,0,4', flags='sf')
Image('view.png')


Out[33]:

In [34]:
grass.run_command('d.mon', stop='png')


Out[34]:
0

In [35]:
# thinned by 125th points
pcloud_sfm_t125 = '_path_to_your_files_/pcloud_sfm_t125.txt'
sfm_t125_density = 'sfm_t125_density'
grass.run_command('g.region', vector=mask_sfm, res=1, flags='pa')

In [36]:
grass.run_command('r.in.xyz', input=pcloud_sfm_t125, output=sfm_t125_density, method='n', separator='comma', overwrite=ow)
grass.run_command('r.null', map=sfm_t125_density, setnull=0)
grass.run_command('r.colors', map=sfm_t125_density, color='viridis', flags='e')

In [37]:
# display with virtual monitor rendered to png
grass.run_command('d.mon', start='png', output='view.png', overwrite=ow)
grass.run_command('d.rast', map=sfm_t125_density)
grass.run_command('d.vect', map=mask_sfm, fill_color='none')
grass.run_command('d.legend', raster=sfm_t125_density, at='5,80,0,4', flags='sf')
Image('view.png')


Out[37]:

In [38]:
grass.run_command('d.mon', stop='png')


Out[38]:
0

In [39]:
# TLS
pcloud_tls_txt = '_path_to_your_files_/pcloud_tls_20mm.txt'
tls_density = 'tls_density'
grass.run_command('g.region', vector=mask_tls, res=1, flags='pa')


Out[39]:
0

In [40]:
grass.run_command('r.in.xyz', input=pcloud_tls_txt, output=tls_density, method='n', separator='comma', overwrite=ow)
grass.run_command('r.null', map=tls_density, setnull=0)
grass.run_command('r.colors', map=tls_density, color='viridis', flags='e')

In [41]:
# display with virtual monitor rendered to png
grass.run_command('d.mon', start='png', output='view.png', overwrite=ow)
grass.run_command('d.rast', map=tls_density)
grass.run_command('d.vect', map=mask_tls, fill_color='none')
grass.run_command('d.legend', raster=tls_density, at='5,60,0,4', flags='sf')
Image('view.png')


Out[41]:

In [42]:
grass.run_command('d.mon', stop='png')


Out[42]:
0

In [43]:
#

In [44]:
###########################################################################
###########################################################################
# 
# Comparison Terrestrial LiDAR / SfM-MVS
# 
###########################################################################
###########################################################################

In [45]:
# Terrestrial LiDAR
# This dataset is available via OpenTopography - 
# exported from FARO Scene as E57, converted to TXT with CloudCompare
# thinned in FARO Scene with 20mm minimum distance between points - 170,141,709 points

In [46]:
# TLS - Terrestrial LiDAR survey
pcloud_tls_txt = '_path_to_your_files_/pcloud_tls_20mm.txt'

In [47]:
# ------------------------------------------------------------------------
# import TLS data as raster (mean val in 10cm), clipped to mask, convert to vector for interpolation
# ------------------------------------------------------------------------
grass.run_command('g.region', vector=mask_tls, res=0.1, flags='pa')
grass.run_command('r.mask', vector=mask_tls, overwrite=ow)
tls_mean = 'tls_rinxyz_mean_10cm'

In [48]:
# import as raster using 'mean'
grass.run_command('r.in.xyz', input=pcloud_tls_txt, output=tls_mean, method='mean', separator='comma', overwrite=ow) 
grass.run_command('r.null', map=tls_mean, setnull=0)

In [49]:
# convert to vector
grass.run_command('r.to.vect', input=tls_mean, output=tls_mean, type='point', flags='zt', overwrite=ow)
# grass.parse_command('v.info', map=tls_mean)
# Number of points:       7,028,118

In [50]:
#

In [51]:
# Sfm-MVS reconstruction of dune field (Agisoft Metashape)
# This dataset is available via OpenTopography - http://dx.doi.org/10.5069/G9DV1H19
# 
# convert LAS to TXT (run this in a terminal):
# cd _path_to_your_files_
# las2txt pcloud_sfm.las --parse xyz --precision 3 3 2 -o pcloud_sfm_full.txt
pcloud_sfm_tlsarea_txt = '_path_to_your_files_/pcloud_sfm_full.txt'

In [52]:
# ------------------------------------------------------------------------
# import SfM data as raster (mean val in 10cm), clipped to mask, convert to vector for interpolation
# ------------------------------------------------------------------------
grass.run_command('g.region', vector=mask_tls, res=0.1, flags='pa')
grass.run_command('r.mask', vector=mask_tls, overwrite=ow)
sfm_mean = 'sfm_rinxyz_mean_10cm'

In [53]:
# import as raster using 'mean'
grass.run_command('r.in.xyz', input=pcloud_sfm_tlsarea_txt, output=sfm_mean, method='mean', separator='comma', overwrite=ow) 
grass.run_command('r.null', map=sfm_mean, setnull=0)

In [54]:
# convert to vector
grass.run_command('r.to.vect', input=sfm_mean, output=sfm_mean, type='point', flags='zt')
# grass.parse_command('v.info', map=sfm_mean)
# Number of points:       8039750

In [55]:
# ------------------------------------------------------------------------
# Interpolate TLS and SfM DEMs, create shaded relief images
# ------------------------------------------------------------------------

In [56]:
# settings for interpolation, shaded reliefs, names for the files
method='bilinear'
step = 0.4
altitude = 30
azimuth = 25
az_txt = '{:>03.0f}'.format(azimuth)

dem_tls_10cm = 'tls_rinxyz_mean_10cm_' + method + '_step_' + str(step)
dem_tls_10cm_shade = dem_tls_10cm + '_shade_' + az_txt + '_' + str(altitude)

dem_sfm_10cm = 'sfm_rinxyz_mean_10cm_' + method + '_step_' + str(step)
dem_sfm_10cm_shade = dem_sfm_10cm + '_shade_' + az_txt + '_' + str(altitude)

diff_sfm_tls_10cm = 'diff_10cm_sfm_tls'
diff_sfm_tls_10cm_shade = diff_sfm_tls_10cm + '_shade_' + az_txt + '_' + str(altitude)

In [57]:
# set region and mask
grass.run_command('g.region', vector=mask_tls, res=0.1, flags='pa')
grass.run_command('r.mask', vector=mask_tls, overwrite=ow)

In [58]:
# interpolate TLS - 10cm
grass.run_command('v.surf.bspline', input=tls_mean, raster_output=dem_tls_10cm, ew_step=step, ns_step=step, method=method, mask='MASK', overwrite=ow)

In [59]:
# interpolate SfM - 10cm
grass.run_command('v.surf.bspline', input=sfm_mean, raster_output=dem_sfm_10cm, ew_step=step, ns_step=step, method=method, mask='MASK', overwrite=ow)

In [60]:
# DEM of Difference - TLS minus SfM
grass.mapcalc('${out} = ${tls} - ${sfm}',
    out=diff_sfm_tls_10cm,
    tls=dem_tls_10cm,
    sfm=dem_sfm_10cm,
    overwrite = True)

In [61]:
# make shaded reliefs and optionally export as tiff
# TLS
grass.run_command('r.relief', input=dem_tls_10cm, output=dem_tls_10cm_shade, azimuth=azimuth, altitude=altitude, overwrite=ow)
# grass.run_command('r.out.gdal', input=dem_tls_shade, output=dem_tls_shade+'.tif', format='GTiff', type='Int16', overwrite=ow)

# SfM
grass.run_command('r.relief', input=dem_sfm_10cm, output=dem_sfm_10cm_shade, azimuth=azimuth, altitude=altitude, overwrite=ow)
# grass.run_command('r.out.gdal', input=dem_sfm_5cm_shade, output=dem_sfm_5cm_shade+'.tif', format='GTiff', type='Int16', overwrite=ow)

# diff 
grass.run_command('r.colors', map=diff_sfm_tls_10cm, color='differences')
grass.run_command('r.relief', input=diff_sfm_tls_10cm, output=diff_sfm_tls_10cm_shade, azimuth=azimuth, altitude=altitude, overwrite=ow)
# grass.run_command('r.out.gdal', input=diff_sfm_tls_shade, output=diff_sfm_tls_shade+'.tif', format='GTiff', type='Int16', overwrite=ow)

In [62]:
# ------------------------------------------------------------------------
# Generate random points, sample DEMs
# ------------------------------------------------------------------------

In [63]:
# import mask to constrain random points 
# (excludes) vegetations and areas without TLS data
shp_tls_sfm = dataDir + 'area_tls_sfm.shp'
mask_tls_sfm = 'mask_tls_sfm'

In [64]:
grass.run_command('v.in.ogr', input=shp_tls_sfm, output=mask_tls_sfm, flags='o', overwrite=ow)

In [65]:
seed = 654321
npoints = 2000
rpoints = 'rpoints'

In [66]:
grass.run_command('v.random', output=rpoints, npoints=npoints, restrict=mask_tls_sfm, seed=seed, overwrite=ow)
grass.run_command('v.db.addtable', map=rpoints, columns="tls double, sfm double, diff double")
grass.run_command('v.what.rast', map=rpoints, raster=dem_tls_10cm, column='tls')
grass.run_command('v.what.rast', map=rpoints, raster=dem_sfm_10cm, column='sfm')
grass.run_command('v.what.rast', map=rpoints, raster=diff_sfm_tls_10cm, column='diff')

In [67]:
# display with virtual monitor rendered to png
grass.run_command('d.mon', start='png', output='view.png', overwrite=ow)
grass.run_command('d.rast', map=dem_tls_10cm_shade)
grass.run_command('d.vect', map=mask_tls, fill_color='none')
grass.run_command('d.vect', map=mask_tls_sfm, fill_color='none', color='red') # mask for random points
grass.run_command('d.vect', map=rpoints)
Image('view.png')


Out[67]:

In [68]:
grass.run_command('d.mon', stop='png')


Out[68]:
0

In [69]:
# ------------------------------------------------------------------------
# Get attr data, make stats and plots
# ------------------------------------------------------------------------

In [70]:
# read vector from GRASS
elevs = gvect.vector_db_select(rpoints)['values']

In [71]:
gvect.vector_db_select(rpoints)['columns']


Out[71]:
['cat', 'tls', 'sfm', 'diff']

In [72]:
# create dataframe, fix data type
# pandas must be newer than 0.22
sfm_tls = pd.DataFrame.from_dict(elevs, orient='index', columns=['cat','tls','sfm','diff'])
sfm_tls['tls'] = sfm_tls['tls'].astype(float)
sfm_tls['sfm'] = sfm_tls['sfm'].astype(float)
sfm_tls['diff'] = sfm_tls['diff'].astype(float)
# sfm_tls.describe()

In [73]:
# calculate stats for each DEM
# could've used pd.describe() but I wanted skew and kurtosis as well
df_stats = pd.DataFrame(columns = ['min','max','mean','median','stddev','skew','kurt','p25','p75'], index=['tls','sfm','diff'])
for dem in ['tls','sfm','diff']:
    rast_max = sfm_tls[dem].max()
    rast_min = sfm_tls[dem].min()
    rast_mean = sfm_tls[dem].mean()
    rast_median = sfm_tls[dem].quantile(q=0.50)
    rast_stddev = sfm_tls[dem].std()
    rast_p25 = sfm_tls[dem].quantile(q=0.25)
    rast_p75 = sfm_tls[dem].quantile(q=0.75)
    rast_skew = sfm_tls[dem].skew()
    rast_kurt = sfm_tls[dem].kurtosis()
    df_stats.loc[dem] = [rast_min,rast_max,rast_mean,rast_median,rast_stddev,rast_skew,rast_kurt,rast_p25,rast_p75]

df_stats
# print df_stats.to_latex(float_format="{:4.2f}".format)
# \begin{tabular}{llllllllll}
# \toprule
# {} &    min &    max &   mean & median & stddev &   skew &   kurt &    p25 &    p75 \\
# \midrule
# tls  &  8.518 & 36.128 & 23.749 & 24.641 &  6.436 & -0.217 & -0.749 & 18.238 & 28.247 \\
# sfm  &  8.684 & 36.104 & 23.849 & 24.796 &  6.410 & -0.231 & -0.755 & 18.393 & 28.327 \\
# diff & -0.623 &  0.226 & -0.100 & -0.103 &  0.129 & -0.259 &  0.119 & -0.196 &  0.007 \\
# \bottomrule
# \end{tabular}


Out[73]:
min max mean median stddev skew kurt p25 p75
tls 8.51795 36.1278 23.7488 24.6413 6.43646 -0.216549 -0.748551 18.2376 28.2471
sfm 8.68434 36.1044 23.849 24.7965 6.41034 -0.230788 -0.755286 18.3929 28.327
diff -0.623076 0.22565 -0.100245 -0.102535 0.129369 -0.259375 0.118626 -0.196481 0.0067143

In [74]:
# how much of the TLS DEM is below the SfM one?
ss.percentileofscore(tuple(sfm_tls['diff']),0)


/usr/local/lib/python2.7/site-packages/scipy/stats/stats.py:1792: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  pct = (np.mean(a_len[idx]) / n) * 100.0
Out[74]:
73.7

In [75]:
# boxplot 
bp = plt.boxplot(sfm_tls['diff'])
# plt.savefig('diff_boxplot.svg')