Vector Operations

A nice introduction to the GRASS vector model can be found here.

Importing and Exporting from/to a Garmin GPS

Importing and exporting GRASS vector data from/to GPS unit can be done in many ways. The commands v.in.garmin and v.in.gpsbabel work well for importing track and waypoint data from a GPS. Note however that v.in.garmin requires points and track data to be stored in decimal degree format. A simpler approach for dealing with both import and export of vector data from GRASS cab be done with a little bit of awk and the program gpstrans. A few common examples are listed below. Command names are in boldface type, input files are red, and the resulting GRASS data are in green.

Gpsbabel

From GRASS

# garmin-specific notes here
http://www.gpsbabel.org/htmldoc-1.3.2/fmt_garmin.html

1. dowload data vis gpsbabel to GPX:
 - import with v.in.gpsbabel
 - convert directly to shp with gps2shp:
   http://gpx2shp.sourceforge.jp/  (source)

Samples from the command line.

# download waypoints from garmin to GPX file
gpsbabel -i garmin -f /dev/ttyS0 -o gpx -F garmin_points.gpx

# download track from garmin to GPX file
gpsbabel -t -i garmin -f /dev/ttyS0 -o gpx -F garmin_track.gpx

# upload to garmin GPS
gpsbabel -i gpx -f garmin_points.gpx -o garmin -F /dev/ttyUSB0

Gpstrans

Import all way points with labels starting with "0":

#download all waypoints from gps
gpstrans -dw > all_wpts_5-3.ascii
#output from gpstrans looks like:

Format: UTM  UTC Offset:  -7.00 hrs  Datum[066]: NAD83
W       001     01-MAY-06 22:58                                 12/31/1989 -7:00:00     10      S       0662251 4039089
W       002     01-MAY-06 23:20                                 12/31/1989 -7:00:00     10      S       0661915 4039282
W       003     01-MAY-06 23:36                                 12/31/1989 -7:00:00     10      S       0662224 4039287
W       004     02-MAY-06 00:06                                 12/31/1989 -7:00:00     10      S       0662002 4039768

#import only those which start with "0" into GRASS
awk '{OFS="|"} $2 ~ /^0.?/{print $9,$10,$2}' all_wpts_5-3.ascii | v.in.ascii out=new_vector columns="x int, y int, id varchar(3)"

Export from points to Garmin GPS:

#export with a "P" prefix to point names

v.out.ascii spring3 dp=0 | \
awk -F"|" '
BEGIN{
OFS="\t"; print "Format: UTM  UTC Offset:  -7.00 hrs  Datum[066]: NAD83"
}
{
print "W\tP"$3"     \t28-MAR-07 20:35                         \t12/31/1989 -7:00:00\t10\tS\t0"$1"\t"$2
}
' > garmin_points.waypts

#upload new waypoints to gps
gpstrans -uw garmin_points.waypts

Export new route to Garmin GPS:

##now make a route file, using points P1 --> P6
v.out.ascii pygmy_route_points \
| awk -F"|" 'BEGIN{OFS="\t"; print "Format: UTM UTC Offset: -7.00 hrs Datum[066]: NAD83\nR\t0\t________________"}
{split($1,easting,"."); split($2,northing,".");
print "W","P"$3,"02-MAY-06 17:19","12/31/1989 -7:00:00","10","S","0"easting[1],northing[1]}' > pygmy_route_points.route

#upload new waypoints to gps
gpstrans -ur pygmy_route_points.route

Raster profile along arbitrary line segments

Overview
A simple experiment illustrating the calculation of elevation profiles along arbitrary line segments using a combination of GRASS and R. Fractal dimension was calculated for each profile by spectral decomposition (FFT), as suggested in (Davis, 2002). Note that the elevation profiles are calculated in a rather non-standard way, as elevation is plotted against 3D distance along the original line segment. This means that resulting transect length calculations are sensitive to the resolution of the raster (or distance between sample points), across a rough surface: i.e. smaller grid spacing will result in a longer effective (3D) transect length.

References

Davis, J.C. Statistics and Data Analysis in Geology John Wiley and Sons, 2002

Sample profile creation I: Vector representation of 2 transects performed by an adventurous soil scientist near McCabe Canyon, Pinnacles National Monument.

Create 2 profiles in GRASS

 # digitize some new lines in an area of interest

v.digit -n map=valley bgcmd="d.rast elev_meters"

v.digit -n map=mtns bgcmd="d.rast elev_meters"

# convert to points

v.to.points -i -v -t in=valley out=valley_pts type=line dmax=50 

v.to.points -i -v -t in=mtns out=mtns_pts type=line dmax=50 

# set resolution to match DEM

g.region -pa res=4 

# extract elevation at each point

v.drape in=valley_pts out=valley_pts_3d type=point rast=elev_4m method=cubic

v.drape in=mtns_pts out=mtns_pts_3d type=point rast=elev_4m method=cubic

# export to text files

v.out.ascii valley_pts_3d > valley_pts_3d.txt

v.out.ascii mtns_pts_3d > mtns_pts_3d.txt

# start R

R

Sample profile creation II: Elevation profiles for the paths taken by the soil scientist here.

Process profile data, and plot in R

# read in the data, note that there is no header, and the columns are pipe delimited

valley <- read.table('valley_pts_3d.txt', header=F, sep="|", col.names=c('easting','northing','elev','f'))

mtns <- read.table('mtns_pts_3d.txt', header=F, sep="|", col.names=c('easting','northing','elev','f'))

# compute the length minus 1 ; i.e. the number of points - 1

# use for the lagged distance calculation

s.valley <- length(valley$easting) - 1

s.mtns <- length(mtns$easting) - 1

# calculate the lagged distance along each component of the 3D vector

dx.valley <- valley$easting[2:(s.valley+1)] - valley$easting[1:s.valley]

dy.valley <- valley$northing[2:(s.valley+1)] - valley$northing[1:s.valley]

dz.valley <- valley$elev[2:(s.valley+1)] - valley$elev[1:s.valley]

dx.mtns <- mtns$easting[2:(s.mtns+1)] - mtns$easting[1:s.mtns]

dy.mtns <- mtns$northing[2:(s.mtns+1)] - mtns$northing[1:s.mtns]

dz.mtns <- mtns$elev[2:(s.mtns+1)] - mtns$elev[1:s.mtns]

# combine the vector components

# 3D distance formula

c.valley <- sqrt( dx.valley^2 + dy.valley^2 + dz.valley^2)

c.mtns <- sqrt( dx.mtns^2 + dy.mtns^2 + dz.mtns^2)

# create the cumulative distance along line:

dist.valley <- cumsum(c.valley[1:s.valley])

dist.mtns <- cumsum(c.mtns[1:s.mtns])

# plot the results

par(mfrow=c(2,1))

plot(dist.valley, valley$elev[1:s.valley], type='l', xlab='Distance from start (m)', ylab='Elevation (m)', main='Profile')

plot(dist.mtns, mtns$elev[1:s.mtns], type='l', xlab='Distance from start (m)', ylab='Elevation (m)', main='Profile')

Sample profile creation III: Spectral decomposition (by FFT) for each profile.

Compute fractal dimension of each transect by FFT

 #compute the fractal dimension via [1]

s.valley <- spectrum(valley$elev[1:s.valley])

s.mtns <- spectrum(mtns$elev[1:s.mtns])

# compute linear model of log(power) ~ log(1/frequency)

# use the slope estimate in [2]

# first define the model formuals: for plotting and for lm()

fm.valley <- log(s.valley$spec) ~ log(1/s.valley$freq)

fm.mtns <- log(s.mtns$spec) ~ log(1/s.mtns$freq)

# next compute the least-squares regression lines via lm()

lm.valley <- lm(fm.valley)

lm.mtns <- lm(fm.mtns)

# extract the slope of each regression line:

slope.valley <- coef(lm.valley)[2]

slope.mtns <- coef(lm.mtns)[2]

# compute Df - the fractal dimension from [3]

Df.valley <- (5 - slope.valley) / 2

Df.mtns <- (5 - slope.mtns) / 2

# plot the figures

plot(fm.valley, main=paste("Valley Profile\nFractal Dimension: ", round(Df.valley,3), sep=''), xlab='Log(1/frequency)', ylab='Log(Power)' )

abline(lm.valley)

plot(fm.mtns, main=paste("Mountain Profile\nFractal Dimension: ", round(Df.mtns,3), sep=''), xlab='Log(1/frequency)',  ylab='Log(Power)')

abline(lm.mtns)

Spatial Clustering of Point Data: Spearfish Example

This example uses the 'Partitioning Around Medoids (PAM)' algorithm (Kaufman and Rousseeuw, 2005) to divide a number of point observation into k clusters, based on their spatial attributes only. An extension of this concept can be applied to any type of geographic data, such as terrain attributes.

For a simple comparison of some of the partitioning-style clustering algorithms in R, see this page of demos. For a more complete listing of clustering approaches in R, see the Cluster Task View.

References

Kaufman, L. & Rousseeuw, P.J. Finding Groups in Data An Introduction to Cluster Analysis Wiley-Interscience, 2005

Export xy coordinates for the bugsites from GRASS See attached file at bottom of page.

# export bugsites

v.out.ascii in=bugsites out=bugsites.xy

Load this text file into an R session A simple map can be made by plotting the xy coordinates.

# read in ascii file, and assign column names

x <- read.table('bugsites.xy', sep="|")

names(x) <- c('easting', 'northing', 'cat')

&nbsp;

# subset original object, return only x,y cols

y <- data.frame(x[,1:2])

row.names(y) <- x$cat

&nbsp;

# simple plot of x,y data

plot(y, pch=3)

Flexclust Example: determining the right number of clusters

Use the stepFlexclust function to determine an optimal number of hard classes 5 clusters looks like a good start.

# load cluster package

library(cluster)

library(flexclust)

&nbsp;

# figure out a good number of clusters use a range of 2 to 10 clusters, with 20 reps each

s <- stepFlexclust(y, K=2:10, nrep=20)

plot(s)

PAM Example: plot of the 5 clusters, centroids, and connectivity

Perform hard classification (clustering) with the PAM algorithm, and plot the results

# 5 clusters in a good number

y.pam <- pam(y, 5, stand=TRUE)

&nbsp;

# plot the clusters by color

plot(y$easting, y$northing, col=y.pam$clustering, main="Bugsites Spatial Clustering, 5 classes", cex=0.5, pch=16, xlab="Easting", ylab="Northing")

&nbsp;

# add the centroids

points(y.pam$medoids, pch=15, col=1:5, cex=1.25)

&nbsp;

# connect the original points to the centroids with line segments:

for(i in 1:5) 

{

segments(x0=y.pam$medoids[i,][[1]], y0=y.pam$medoids[i,][[2]], x1=y$easting[y$cluster == i], y1=y$northing[y$cluster ==i], col=i, lty=3)

}

Prepare the data for export to text, and save the clustered data

# add the cluster number to the original dataframe

y$cluster <- y.pam$clustering

y$orig_cat <- as.numeric(row.names(y))

&nbsp;

# save as a text file and quit

write.table(y, file='bugsites.clust', row.names=FALSE)

Load clustered data into GRASS as a new set of points called 'bclust' For each cluster, extract those points, and compute a convex hull.

# load clustered points into GRASS

v.in.ascii in=bugsites.clust out=bclust fs=" " columns='x double, y double, cluster integer, orig_cat integer' skip=1

&nbsp;

# zoom to the full extent of the Spearfish dataset

g.region rast=elevation.dem

&nbsp;

# there are 5 clusters: show them all, and compute convex hulls

for x in `seq 1 5` 

do v.extract --o in=bclust where="cluster=$x" out=bclust_$x 

v.hull --o in=bclust_$x out=bclust_hull_$x

d.vect bclust_hull_$x type=boundary fcol=none width=2 col=white

d.vect bclust icon=basic/box fcol=black col=black size=6

done

Traveling Salesman Approach to Visiting Data-loggers

Optimal Route in 3D

Premise:
Several data-loggers in steep terrain; how can we visit all of them in the most efficient manner such that we traverse the least cumulative slope gradient? In a sense this is a traveling salesman style problem, with travel costs defined by slope gradient. The v.net.salesman module in GRASS can be used to compute the optimal route between data-loggers, given: 1) a network of lines connecting all points to be visited and 2) transferal of slope cost information to the network. An example approach is listed below. Note that due to a bug in v.net.salesman, we need to use one extra step-- details below.

Approach:

generate a delaunay network between points
densify the line segments along the network (v.split)
convert the network to a 3D, by sampling a DEM along vertices (v.drape)
save the slope of each line segment on the network to its attribute table (v.to.db)
use the squared slope as the travel cost
re-connect the new attribute table from the 3D map, to the original 2D map (v.db.connect)
connect data-logger points to delaunay network (v.net)
solve the traveling salesman problem (v.net.salesman)

Room for Improvement:

generate a cost surface using a combination of slope map and known walkable trails / roads
add support for lines to v.what.rast -- this would allow us to sample a more realistic cost surface, generated with r.cost
fix v.net.salesman so that it can use 3D input maps

Optimal Datalogger Route

Implementation in GRASS:

# generate simple network connecting all points

v.delaunay --o -l in=temperature_sites_maybe out=tri

# densify line segments

v.split --o in=tri out=tri_2 length=5

# add categories

v.category --o in=tri_2 out=tri_3 option=add type=line

# convert to 3d, using our elevation model

v.drape --o in=tri_3 out=tri_4 rast=rtk_elev_var_smooth@rtk

# add a table

v.db.addtable map=tri_4 columns="cat integer, cost double"

# use the slope between line segments as the cost

v.to.db map=tri_4 type=line layer=1 option=slope columns=cost

# make all slopes positive

echo "UPDATE tri_4 set cost = cost * cost" | db.execute

# connect table from 3D map to original 2D map:

v.db.connect map=tri_3 table=tri_4

# make a copy, which brings over a copy of the linked table

g.copy --o vect=tri_3,triangulation

# clean-up

g.remove vect=tri,tri_2,tri_3,tri_4

# make a network

v.net --o in=triangulation points=temperature_sites_maybe out=net --o operation=connect thresh=10

# check network-- has categories?

v.category net op=report

# solve traveling salesman problem

# does not work with 3D maps (ticket # 406)

v.net.salesman --o in=net out=salesman acol=cost ccats=1-33

Traveling Salesman Approach to Visiting Data-loggers II

Premise:
Updated version of a previous attempt to use a traveling-salesman approach to visiting data-loggers in steep terrain. This time we use a vector/raster approach (v.rast.stats) to compute the total "cost" incurred along each segment of a network connecting all of our data-loggers. Although this approach requires more computation time (v.rast.stats is slow), we can generate a more detailed cost surface to model movement. In this case we are using a cumulative-traversed-slope derived cost, with a 90% reduction along a local road.

New Approach:

generate a cumulative cost surface, originating at data-logger locations, based on slope angle (r.slope.aspect, r.cost)
generate a network of lines connecting all data-loggers, along with a local road (v.delaunay, v.digit, v.patch)
modify the cost surface along the local road, such that cumulative cost is reduced 90% within 5m of the road (r.buffer, r.mapcalc)
densify the line segments along the network (v.split)
compute the total cost along each line segment of our network (v.rast.stats)
connect data-logger points to delaunay network (v.net)
solve the traveling salesman problem (v.net.salesman)

Optimal Route version 2: 10-meter contours in grey, network in blue, optimal cycle between data-loggers (yellow) in red.

Setup the cost surface and network

# generate cost surface based on slope

r.slope.aspect elev=rtk_elev_var_smooth@rtk slope=slope format=percent

r.cost --o -k in=slope start_points=temperature_sites_maybe output=slope_cost

# make a vector of roads within the watershed:

v.digit -n map=roads  bgcmd="d.rast naip2005.green"

# use this to modify the cost surface

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

# generate buffer 5 meters from road

r.buffer --o in=roads out=roads_buff distances=5 units=meters

# reduce slope cost to 1% of its value along the buffered road

r.mapcalc "new_cost = if(isnull(roads_buff), slope_cost, (0.01) * slope_cost)"

# generate simple network connecting all points

v.delaunay --o -l in=temperature_sites_maybe out=tri

# combine road network, with delaunay network

v.patch in=tri,roads out=tri_1

v.clean --o in=tri_1 out=tri_1_clean tool=break thresh=1

# densify line segments

v.split --o in=tri_1_clean out=tri_2 length=20

# add categories

v.category --o in=tri_2 out=tri_3 option=add type=line

# add table

v.db.addtable map=tri_3 columns="cat integer"

Compute the total cost along each line segment, solve traveling salesman problem

# better approach to assigning cost, based on cost surface:

# really slow, as it uses vect->rast conversion to compute statistics

# about 3.5 minutes for 588 line features =  features / minute

# set the region to 2m, so that rasterized lines are wider

g.region res=2

v.rast.stats vect=tri_3 layer=1 rast=new_cost colprefix=s

# reset region

g.region res=1

# make network

v.net --o in=tri_3 points=temperature_sites_maybe out=net --o operation=connect thresh=10

# check network-- has categories?

v.category net op=report

# solve traveling salesman problem -- does not work with 3D maps

# use regular path-distance cost:

v.net.salesman --o in=net out=simple_salesman ccats=1-33

# use the total cost surface values / line segment

v.net.salesman --o in=net out=salesman acol=s_sum ccats=1-33

Optimal Route with cost surface: optimal route in black, cost surface colors: blue = low, red = high.

Working with transects

Figure 1: screen plotting.

Figure 2: new vector creation.

Figure 3: points along line.

Transects are often used to better understand local variation in soil properties. Planning and mapping a transect in the office can make time spent in the field more efficient. Sometimes we return from the field with a transect starting point and bearing (CCW from north) of travel. Accurate digitizing of this information can often be time consuming when multiple transects are involved.

The command d.bearing can facilitate rapid planning, visualization, and digitization of transects given the parameters: origin=x,y, theta=bearing, radius=length_of_transect. The following examples illustrate some common methods of working with transect data. Command names are in boldface and comments are in grey. The source code for d.bearing is attached at the bottom of this page. Unarchive in the grass6/display/ directory of the GRASS61-CVS source code directory.

Example 1: Plot a hypothetical transect; bearing 230 degrees, 1000m length.

Setup the display:

#plot the transect origin point:
d.vect pedons cats=91 icon=basic/box
#make a grid for illustative purposes, with origin at the point above:
d.grid size=100 -t -b origin=665736,4039451
#add a quick scale bar with the mouse:
d.barscale -m

Plot the transect to the screen. (Figure 1)

d.bearing origin=665736,4039451 theta=230 radius=1000 col=blue

Example 2: Make a permanent vector 't1' of the same transect.

Replot the transect, this time saving the transect line to a new vector called 't1'

d.bearing origin=665736,4039451 theta=230 radius=1000 col=blue output=t1

Setup the display again, this time plotting the new transect vector 't1'. (Figure 2)

d.erase
d.grid size=100 -t -b origin=665736,4039451
d.vect pedons cats=91 icon=basic/box
#display the new vector 't1' along with its direction
d.vect t1 col=blue disp=shape,dir

Example 3: Create a new vector of points every 100m along the transect 't1'.

Assign unique IDs to the vertices of transect vector 't1'

#setup a database connection for vector 't1' to MySQL
v.db.connect map=t1 driver=mysql database="host=localhost,dbname=test" table=t1 key=cat
#create a new table in MySQL to store information about vector 't1'
echo "create table t1 (cat int not null , primary key (cat) )" | db.execute
#upload a unique ID of each vertex to the new table for 't1'
v.to.db map=t1 option=cat

Extract a set of points to 't1_points' every 100m along transect vector 't1'.

#extract a new point vector along 't1' interpolating a new point every 100m, saving to 't1_points'
v.to.points -i in=t1 out=t1_points dmax=110
#the new vector 't1_points' now has two tables associated with it:

t1_points_1: layer 1

t1_points_2: layer 2

+------+
| cat  |
+------+
|    1 |
|    2 |
+------+

+------+------+------------------+
| cat  | lcat | along            |
+------+------+------------------+
|    1 |    1 |                0 |
|    2 |    1 | 100.000000000004 |
|    3 |    1 | 200.000000000007 |
|    4 |    1 | 300.000000000011 |
|    5 |    1 | 400.000000000015 |
|    6 |    1 | 500.000000000019 |
|    7 |    1 | 600.000000000022 |
|    8 |    1 | 700.000000000026 |
|    9 |    1 |  800.00000000003 |
|   10 |    1 | 900.000000000034 |
|   11 |    1 | 1000.00000000004 |
+------+------+------------------+

display our new point vector 't1_points', labeling the new ID numbers stored in layer 2. (Figure 3)

d.vect t1_points icon=basic/point size=10 fcol=red col=red disp=shape,cat llayer=2