Elevation Survey by RTK

RTK base station setup

Draft Elevation Model

Vertical Precision: Yellow points represent vertical precisions less than 10cm, red points are for precisions greater than 10cm. 10 meter contour interval.

Regions with low point coverage: in yellow. 10 meter contour interval.

Draft SFREC elevation model: viewed from the north-west

Estimates of Local Sky Obstruction (Regions of Problematic GPS Reception)

Some possible approaches

r.angle.sh (see image below and attached script)
r.mapcalc "$GIS_OPT_output = atan( ($GIS_OPT_elev - $z) / sqrt( (x() - $x)^2 + (y() - $y)^2) )"

Angle Map: Vertical angle (degrees) from starting point (yellow point) to each cell.
Compute percent of visible terrain that is above 10 degrees vertical elevation at each cell, using r.los.
# setup:
#
# 30 m resolution takes about 4.5 minutes on my laptop
# g.region res=30 -ap
# 10 m resolution takes about 19 minutes on my laptop
# g.region res=10 -ap

# dump x,y,z values from DEM
r.out.xyz elev10m > ~/temp/xyz

# setup coordinate file
echo -n "" > stats.xyp
g.remove --q rast=MASK

# loop through each cell
for cell in `cat ~/temp/xyz`
do

los_cmd=`echo $cell | awk -F"|" '{print "r.los --q --o in=elev10m out=los coordinate="$1","$2}'`

angle_cmd=`echo $cell | awk -F"|" '{print "r.mapcalc \"ang = if( los >= 100, 1, null()) \""}'`

eval $los_cmd
eval $angle_cmd

p=`r.stats -p -N ang | awk '{print $2}' | sed 's/\%//g'`
echo $cell | awk -v p=$p -F"|" '{ if(p != "") {print $1"|"$2"|"p}}' >> stats.xyp

done

# load results when done:
# r.in.xyz in=stats.xyp out=angp --o

Percent of visible terrain obscuring the sky above 10 degrees

Preliminary Watershed Walk

Track Log Data (white flow lines are from original 10 meter USGS DEM)
Interpolated elevation map from Garmin 120 track logs

Location of Garmin track points

Establishing the Optimal Number of Survey Points

First Pass: Terrain Skeleton

Terrain Skeleton

Terrain Skeleton 3D

Initial Estimation

50 meter interval sampling of contours
10 meter interval sampling of terrain skeleton
15 meter thinning (snapping) radius

Lewis 1 RTK Sampling Plan 832 points

# generate a slightly larger watershed, so that we get the edges

v.buffer in=w_smooth out=w_buff buffer=50 --o

# setup some kind of sampling intervals

contour_sampling_interval=50

skeleton_sampling_interval=15

# create a set of points along contours

v.to.points -i -t in=contours out=cp dmax=$contour_sampling_interval --o

v.category in=cp out=cp_cat option=del --o

v.category in=cp_cat out=cp_cat_1 option=add --o

# select only those within the watershed:

v.select -t ainput=cp_cat_1 atype=point binput=w_buff btype=area operator=overlap out=contour_pts --o

# create set of points from terrian skeleton

# note that the output from r.flow involves multiple line segments that converge

# and thus is not readily converted to a regular interval of points

# combine flow and upflow

v.patch in=flow,upflow out=skel

# convert to points on a regular interval

# use a small number to preserve the complexity-- 

# we will thin the number of points down later

v.to.points -i -t in=skel out=sp dmax=10 --o

# these points need a category first

# 

v.category in=sp out=sp_cat option=del --o

v.category in=sp_cat out=sp_cat_1 option=add --o

# select only points within watershed 

v.select -t ainput=sp_cat_1 atype=point binput=w_buff btype=area operator=overlap out=skel_pts --o

# give each sampling point a unique ID and attr designating the observation type:

v.db.addtable skel_pts

v.db.addtable contour_pts

v.db.addcol skel_pts column='pt_type varchar(10)'

v.db.addcol contour_pts column='pt_type varchar(10)'

echo "update contour_pts set pt_type = 'contour'" | db.execute

echo "update skel_pts set pt_type = 'skel'" | db.execute

# patch contour and skeleton sampling points

v.patch -e in=contour_pts,skel_pts out=rp

v.build rp

# thin points out a bit

v.clean in=rp out=rtk_pts tool=snap thres=$skeleton_sampling_interval --o

# cleanup temp points:

g.remove vect=cp,cp_cat,cp_cat_1,sp,sp_cat,sp_cat_1,rp

# simple map

d.erase

d.vect contours col=grey

d.vect skel col=grey

d.vect w_buff width=2 type=boundary

# sample points

# d.vect contour_pts icon=basic/box fcol=yellow

# d.vect skel_pts icon=basic/box fcol=red

# combined sample points:

d.vect rtk_pts icon=basic/box fcol=yellow

# d.barscale -m

# d.out.file intial_plan --o

# how many points are we talking about ?

eval `v.info rtk_pts -t | grep points`

echo $points

# 832

Regions of like SD(elevation)

Compute SD of elevation on 50x50 meter grid

# prepage SD map

r.neighbors in=e10 out=e10_sd method=stddev size=5 --o

Cluster on (standardized) coordinates and SD of elevation

library(cluster)

library(spgrass6)

x <- readRAST6('e10_sd')

n_classes <- 5

x.clara <- clara(data.frame(coordinates(x),  slot(x, 'data')), k=n_classes, stand=TRUE)

x$clust <- factor(x.clara$clustering)

spplot(x, zcol='clust', col.regions=terrain.colors(n_classes))

# pdf(file='elev_sd_clust.pdf', height=5, width=5)

boxplot(x$e10_sd ~ factor(x$clust, levels=order(tapply(x$e10_sd, x$clust, median)) ), ylab='Standard Deviation (meters)', xlab='Class', main='Standard Deviation of Elevation on 50x50 Grid', cex=0.5 )

# dev.off()

x$clust <- x.clara$clustering

writeRAST6(x, zcol='clust', vname='t_clust')

Finish up in GRASS

# convert to int map

r.mapcalc "t_clust = int(t_clust)"

# vectorize

# r.to.vect t_clust out=t_clust feature=area --o

# simple map

d.rast.leg t_clust

d.vect contours

d.vect upflow col=white

d.vect flow col=white

# save

d.out.file out=sd_clust_and_flow

Regions of like SD(elevation)

Elevation standard deviation classes

Selection by Minimizing MAE (mean absolute error)

5 meter interval sampling of contours
5 meter interval sampling of terrain skeleton
[5 10 15 20 25 30 40 50 60] meter thinning (snapping) radius

MAE vs. Thinning radius: terrain skeleton approach in red

# # prep work:

g.region res=10

v.to.rast in=w_buff out=w_buff_mask use=val val=1 --o

# invoked like this

# time sh eval_pts_spacing.sh 2>/dev/null

# approx 7 mins

# start logfie

outfile='stats.out'

echo '' > $outfile

# the min distance between new points when converting lines to points

# make this a small number for lots of points

dmax=5

# sequence of thinning thresholds

sequence="5 10 15 20 25 30 40 50 60"

#       loop over final cleaning snapping threhold

for j in $sequence

        do

        # feedback

        echo "thinning at :  $j meters"

        # setup some kind of sampling intervals

        skeleton_sampling_interval=$j

        # create a set of points along contours

        v.to.points --q -i -t in=contours out=cp_00001 dmax=$dmax --o

        v.category --q in=cp_00001 out=cp_cat_00001 option=del --o

        v.category --q  in=cp_cat_00001 out=cp_cat_1_00001 option=add --o

        # select only those within the watershed:

        v.select --q --o -t ainput=cp_cat_1_00001 atype=point binput=w_buff btype=area operator=overlap out=contour_pts_00001

        # convert to points on a regular interval

        # use a small number to preserve the complexity-- 

        # we will thin the number of points down later

        v.to.points --q  -i -t in=skel out=sp_00001 dmax=$dmax --o

        # these points need a category first

        # 

        v.category  --q in=sp_00001 out=sp_cat_00001 option=del --o

        v.category  --q in=sp_cat_00001 out=sp_cat_1_00001 option=add --o

        # only points within watershed

        v.select -t --o --q ainput=sp_cat_1_00001 atype=point binput=w_buff btype=area operator=overlap out=skel_pts_00001

        # give each sampling point a unique ID and attr designating the observation type:

        v.db.addtable  --q --o skel_pts_00001 2>/dev/null

        v.db.addtable  --q --o contour_pts_00001 2>/dev/null

        v.db.addcol  --q --o skel_pts_00001 column='elev double'

        v.db.addcol  --q --o contour_pts_00001 column='elev double'

        # patch contour and skeleton sampling points

        v.patch  --o --q -e in=contour_pts_00001,skel_pts_00001 out=rp_00001

        v.build  --q rp_00001

        # thin points out a bit

        v.clean  --o --q in=rp_00001 out=rtk_pts_00001 tool=snap thres=$skeleton_sampling_interval

        # extract number of points:

        # save to $points

        eval `v.info rtk_pts_00001 -t | grep points`

        echo "$points points along terrain skeleton"

        #extract elevation at each sampling points

        v.what.rast --q vect=rtk_pts_00001 raster=e10 column=elev

        # interpolate elevation for current region

        # at the same time capture the number of points used in the interpolation

        interp_points=`v.surf.rst in=rtk_pts_00001 zcol=elev elev=e_interp_00001 maskmap=w_buff_mask --o | grep 'The number of points from vector map is' | awk -F'is ' '{print $2}'`

        # generate n random points, where n= number of points used in terrain skeleton-based interpolation

        # note that these are generated along the grid used for sampling

        # only generate points within watershed mask

        r.random --o in=e10 n=$interp_points vector=rand_00001 cover=w_buff_mask

        # interpolate elevation for current region      (random points) 

        v.surf.rst in=rand_00001 zcol=value elev=rand_interp_00001 maskmap=w_buff_mask --o

        # compute MAE

        r.mapcalc "e_diff_00001 = abs(e10 - e_interp_00001)"

        mae=`r.univar e_diff_00001 | grep 'mean:' | awk -F': ' '{print $2}'`

        # compute MAE

        r.mapcalc "rand_diff_00001 = abs(e10 - rand_interp_00001)"

        rand_mae=`r.univar rand_diff_00001 | grep 'mean:' | awk -F': ' '{print $2}'`

        # save stats

        echo "$j,$points,$interp_points,$mae,$rand_mae" >> $outfile

done

	thinning	points	interp_points	mae	rand_mae
1	5	3375	3318	0.59	0.12
2	10	1730	1693	0.60	0.28
3	15	1058	1040	0.60	0.41
4	20	698	691	0.66	0.57
5	25	501	495	0.67	0.76
6	30	381	379	0.76	1.03
7	40	231	227	0.91	1.89
8	50	152	150	1.17	2.15
9	60	117	117	1.35	2.61