Generation of Vegetation Sampling Areas
Premise:
We need some randomly selected areas for performing a vegetation survey. Two possible approaches are given: random selection of grid cells from a 7x7 grid superimposed over the region of interest, and random generation of 100x100 meter blocks which meet certain criteria on how close they can be from each other and the edges of the region. The region and grid files are attached at the bottom of this page.
Generate Region Data
# this ensures that we always get a cell that is completely in the region
v.mkgrid map=grid position=region grid=7,7 --o
# simple polygon defining current region
v.in.region out=reg
v.mkgrid map=grid position=region grid=7,7 --o
# simple polygon defining current region
v.in.region out=reg
Compute Random Grid Cells
## load libs
library(maptools)
library(spgrass6)
library(spatstat)
## approach 1: select random grid cell
## read in data
x <- readVECT6('grid')
## note the comma in there, to select a named row and all columns
grid10 <- x[sample(length(slot(x, "polygons")), 10), ]
## write back to GRASS (or where ever)
writeVECT6(grid10, 'grid10')
library(maptools)
library(spgrass6)
library(spatstat)
## approach 1: select random grid cell
## read in data
x <- readVECT6('grid')
## note the comma in there, to select a named row and all columns
grid10 <- x[sample(length(slot(x, "polygons")), 10), ]
## write back to GRASS (or where ever)
writeVECT6(grid10, 'grid10')
Random Grid Cell Example
Generate Randomly Placed Blocks
## load libs
library(maptools)
library(spgrass6)
library(spatstat)
## read in the rectangle which defines in the region of interest
reg <- readVECT6('reg')
## load custom function (see attached file below...)
source('functions.R')
## use custom function (see attached file below...)
## spatstat style sampling: with rSSI()
SP <- rnd_blocks(region=reg, n.blocks=10, block.edge.length=100, min.rad.sep=2, min.rad.edge=1)
## save result back to GRASS
writeVECT6(SP, 'blocks')
library(maptools)
library(spgrass6)
library(spatstat)
## read in the rectangle which defines in the region of interest
reg <- readVECT6('reg')
## load custom function (see attached file below...)
source('functions.R')
## use custom function (see attached file below...)
## spatstat style sampling: with rSSI()
SP <- rnd_blocks(region=reg, n.blocks=10, block.edge.length=100, min.rad.sep=2, min.rad.edge=1)
## save result back to GRASS
writeVECT6(SP, 'blocks')
Random Block Example
Additional Example: Compute 3 blocks per "half" of the watershed.
Random Block Example 2
Conversion of Sampling Blocks to GPS Waypoints
Conversion of Sample Plot Boundaries to Points
Example
# extract vertices and centroid
# interpolate vertices resulting in a vertex every 50m
v.to.points -v -i in=blocks_1 out=b1p dmax=75 --o
v.to.points -v -i in=blocks_2 out=b2p dmax=75 --o
# patch them together:
v.patch in=b1p,b2p out=vp --o
v.build vp
# remove duplicate points
# new file does not have an attribute table
v.out.ascii vp | sort -g | uniq | v.in.ascii -t out=veg_points
# add attribute table:
v.db.addtable veg_points
# cleanup
g.remove vect=b1p,b2p,vp
# export to shapefile:
v.out.ogr -e in=veg_points type=point dsn=veg_points.shp
# convert to GPX [LL WGS84]
# translate the 'cat' column to 'name'
ogr2ogr -t_srs '+proj=latlong +datum=WGS84' \
-f GPX -dsco GPX_USE_EXTENSIONS=YES \
veg_points.gpx veg_points.shp \
-sql "SELECT cat AS name FROM veg_points"
# upload the wp_gps.gpx file to the GPS unit
gpsbabel -i gpx -f veg_points.gpx -o garmin -F /dev/ttyS0
# simple PDF map:
# results in a rather large PDF
# looks like Cairo is used by default
d.font.freetype font=VeraBd
d.rast doqq
d.vect veg_points disp=shape,cat icon=basic/circle size=9 lsize=9 fcol=yellow col=black lcol=white
d.out.file format=pdf out=veg_survey/simple_map paper=us-letter
# interpolate vertices resulting in a vertex every 50m
v.to.points -v -i in=blocks_1 out=b1p dmax=75 --o
v.to.points -v -i in=blocks_2 out=b2p dmax=75 --o
# patch them together:
v.patch in=b1p,b2p out=vp --o
v.build vp
# remove duplicate points
# new file does not have an attribute table
v.out.ascii vp | sort -g | uniq | v.in.ascii -t out=veg_points
# add attribute table:
v.db.addtable veg_points
# cleanup
g.remove vect=b1p,b2p,vp
# export to shapefile:
v.out.ogr -e in=veg_points type=point dsn=veg_points.shp
# convert to GPX [LL WGS84]
# translate the 'cat' column to 'name'
ogr2ogr -t_srs '+proj=latlong +datum=WGS84' \
-f GPX -dsco GPX_USE_EXTENSIONS=YES \
veg_points.gpx veg_points.shp \
-sql "SELECT cat AS name FROM veg_points"
# upload the wp_gps.gpx file to the GPS unit
gpsbabel -i gpx -f veg_points.gpx -o garmin -F /dev/ttyS0
# simple PDF map:
# results in a rather large PDF
# looks like Cairo is used by default
d.font.freetype font=VeraBd
d.rast doqq
d.vect veg_points disp=shape,cat icon=basic/circle size=9 lsize=9 fcol=yellow col=black lcol=white
d.out.file format=pdf out=veg_survey/simple_map paper=us-letter
Processing Transect Data
Input File Format (see attached CSV file)
start stop species transect block transect.length notes
0 27 open t01 b1 344
68 93 open t01 b1 344
137 148 open t01 b1 344
160 177 open t01 b1 344
209 219 open t01 b1 344
237 318 open t01 b1 344
332 344 open t01 b1 344
27 68 blue oak t01 b1 344
93 137 blue oak t01 b1 344
148 160 blue oak t01 b1 344
177 209 blue oak t01 b1 344
219 237 blue oak t01 b1 344
318 332 blue oak t01 b1 344
[...]
Setup
## load libs
library(lattice)
## read in the data
## note that this contains multiple blank lines
x.data <- read.csv('sfrec_canopy_coverage.csv')
## remove a mistake: there is overlap between black oak and open cover types
x.data <- x.data[-125,]
library(lattice)
## read in the data
## note that this contains multiple blank lines
x.data <- read.csv('sfrec_canopy_coverage.csv')
## remove a mistake: there is overlap between black oak and open cover types
x.data <- x.data[-125,]
Example visualization (see attached PDF)
dotplot(species ~ start + stop | transect, data=x.data,
layout=c(3,6), as.table=TRUE, subscripts=TRUE,
xlab='Transect Distance (ft)', ylab='Cover Type',
key=list(columns=2, text=list(c('Measured', 'Estimated')), lines=list(lty=c(1,1), col=c(1,2))),
panel=function(x, y, subscripts, groups, ...)
{
## plot the points, and setup the graph
panel.dotplot(x, y, subscripts=subscripts, groups=groups, pch=NA, ...)
## convert factor level to number
y_num <- as.numeric(y)
## create a table of the distance, factor level, and start/stop flag
d <- data.frame(x=x, y=y_num, groups=groups[subscripts])
## get the start points
xy_i <- subset(d, select=c(x,y), subset=groups=='start')
## get the stop points
xy_f <- subset(d, select=c(x,y), subset=groups=='stop')
##
## make a line color, based on the notes field
## first level = black
## second level = red (estimated)
lcol=as.numeric(x.data$notes[subscripts])
## plot the lines
panel.segments(xy_i$x, xy_i$y, xy_f$x, xy_f$y, lwd=2, col=lcol)
}
)
layout=c(3,6), as.table=TRUE, subscripts=TRUE,
xlab='Transect Distance (ft)', ylab='Cover Type',
key=list(columns=2, text=list(c('Measured', 'Estimated')), lines=list(lty=c(1,1), col=c(1,2))),
panel=function(x, y, subscripts, groups, ...)
{
## plot the points, and setup the graph
panel.dotplot(x, y, subscripts=subscripts, groups=groups, pch=NA, ...)
## convert factor level to number
y_num <- as.numeric(y)
## create a table of the distance, factor level, and start/stop flag
d <- data.frame(x=x, y=y_num, groups=groups[subscripts])
## get the start points
xy_i <- subset(d, select=c(x,y), subset=groups=='start')
## get the stop points
xy_f <- subset(d, select=c(x,y), subset=groups=='stop')
##
## make a line color, based on the notes field
## first level = black
## second level = red (estimated)
lcol=as.numeric(x.data$notes[subscripts])
## plot the lines
panel.segments(xy_i$x, xy_i$y, xy_f$x, xy_f$y, lwd=2, col=lcol)
}
)
Convert to matrix representation
d <- list()
for(j in levels(x.data$transect))
{
## work with a subset of the data:
y <- subset(x.data, subset=transect == j)
## init a matrix to hold the transect data: wide format
## fill with 0's
z <- matrix(0, ncol=length(levels(y$species)), nrow=max(y$stop))
## for each level of species, populate the corresponding cells of the matrix
for(i in 1:nrow(y))
{
## increment the start by one (shrinking the number of cells by 1)
y_start <- y$start[i] + 1
y_stop <- y$stop[i]
y_col <- as.numeric(y$species[i])
## encode the canopy type
## using powers of 2
z[y_start:y_stop, y_col] <- 2^y_col
}
eval(parse(text=paste('d$', j, ' <- z', sep='')))
}
for(j in levels(x.data$transect))
{
## work with a subset of the data:
y <- subset(x.data, subset=transect == j)
## init a matrix to hold the transect data: wide format
## fill with 0's
z <- matrix(0, ncol=length(levels(y$species)), nrow=max(y$stop))
## for each level of species, populate the corresponding cells of the matrix
for(i in 1:nrow(y))
{
## increment the start by one (shrinking the number of cells by 1)
y_start <- y$start[i] + 1
y_stop <- y$stop[i]
y_col <- as.numeric(y$species[i])
## encode the canopy type
## using powers of 2
z[y_start:y_stop, y_col] <- 2^y_col
}
eval(parse(text=paste('d$', j, ' <- z', sep='')))
}
Matrix representation and canopy overlap simplification
## generate an example
cbind(d$t10[300:320,], rowSums(d$t10[300:320,]))
cbind(d$t10[300:320,], rowSums(d$t10[300:320,]))
*
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 8 0 8
[2,] 0 0 8 0 8
[3,] 0 0 8 0 8
[4,] 2 0 8 0 10
[5,] 2 0 8 0 10
[6,] 2 0 8 0 10
[7,] 2 0 8 0 10
[8,] 2 0 8 0 10
[9,] 2 0 8 0 10
[10,] 2 0 8 0 10
[11,] 0 0 8 0 8
[12,] 0 0 8 0 8
[13,] 0 0 8 0 8
[14,] 0 0 8 0 8
[15,] 0 0 8 0 8
[16,] 0 0 8 0 8
[17,] 0 0 0 16 16
[18,] 0 0 0 16 16
[19,] 0 0 0 16 16
[20,] 0 0 0 16 16
[21,] 0 0 0 16 16
Convert codes
## generate the canopy cover combination table
## using powers of 2
## note that we are leaving out canopy type '5' (open), as there should be no overlap
g <- t(combn(2^(1:3), 2))
g.lookup <- data.frame(apply(g, 2, function(i) levels(x.data$species)[logb(i, base=2)]), code=rowSums(g))
g.lookup.overlap <- data.frame( canopy_type=paste(g.lookup$X1, g.lookup$X2, sep=' / '), code=g.lookup$code)
## now the lookup table to non-overlapping regions
g.lookup.no_overlap <- data.frame(canopy_type=levels(x.data$species), code=2^(1:length(levels(x.data$species))))
## combine
g.lookup.final <- rbind(g.lookup.no_overlap, g.lookup.overlap)
## for each transect compute the linear totals for each canopy type, including overlap
t_sums <- lapply(d, function(i) table(rowSums(i)) )
## re-create the table with the correct canopy type for each transect
t_sums <- lapply(t_sums, function(i) data.frame(canopy=g.lookup.final$canopy_type[match(names(i), g.lookup.final$code)], t_part=as.vector(i)) )
## using powers of 2
## note that we are leaving out canopy type '5' (open), as there should be no overlap
g <- t(combn(2^(1:3), 2))
g.lookup <- data.frame(apply(g, 2, function(i) levels(x.data$species)[logb(i, base=2)]), code=rowSums(g))
g.lookup.overlap <- data.frame( canopy_type=paste(g.lookup$X1, g.lookup$X2, sep=' / '), code=g.lookup$code)
## now the lookup table to non-overlapping regions
g.lookup.no_overlap <- data.frame(canopy_type=levels(x.data$species), code=2^(1:length(levels(x.data$species))))
## combine
g.lookup.final <- rbind(g.lookup.no_overlap, g.lookup.overlap)
## for each transect compute the linear totals for each canopy type, including overlap
t_sums <- lapply(d, function(i) table(rowSums(i)) )
## re-create the table with the correct canopy type for each transect
t_sums <- lapply(t_sums, function(i) data.frame(canopy=g.lookup.final$canopy_type[match(names(i), g.lookup.final$code)], t_part=as.vector(i)) )
Re-attach transect and block id
## convert to dataframe by "row-binding"
t_sums.df <- do.call('rbind', t_sums)
## re-add the transect id
t_sums.df$transect <- substr(row.names(t_sums.df), 1, 3)
## make a lookup table containing transect -> block relationship
t_b.lookup <- unique(subset(x.data, select=c(transect, block, transect.length)))
## join block data
t_b_sums.df <- merge(x=t_sums.df, y=t_b.lookup)
t_sums.df <- do.call('rbind', t_sums)
## re-add the transect id
t_sums.df$transect <- substr(row.names(t_sums.df), 1, 3)
## make a lookup table containing transect -> block relationship
t_b.lookup <- unique(subset(x.data, select=c(transect, block, transect.length)))
## join block data
t_b_sums.df <- merge(x=t_sums.df, y=t_b.lookup)
Compute percent cover by transect
pct_cover_by_transect <- sweep(tapply(t_b_sums.df$t_part, list(t_b_sums.df$canopy, t_b_sums.df$transect), sum, na.rm=TRUE), 2, t_b.lookup$transect.length, '/') * 100
write.csv(pct_cover_by_transect, na='', file='pct_cover_by_transect.csv')
print(pct_cover_by_transect, digits=1)
write.csv(pct_cover_by_transect, na='', file='pct_cover_by_transect.csv')
print(pct_cover_by_transect, digits=1)
Compute percent cover by block
pct_cover_by_block <- sweep(tapply(t_b_sums.df$t_part, list(t_b_sums.df$canopy, t_b_sums.df$block), sum, na.rm=TRUE), 2, tapply(t_b.lookup$transect.length, t_b.lookup$block, sum), '/') * 100
write.csv(pct_cover_by_block, na='', file='pct_cover_by_block.csv')
print(pct_cover_by_block, digits=1)
write.csv(pct_cover_by_block, na='', file='pct_cover_by_block.csv')
print(pct_cover_by_block, digits=1)
Compute percent cover across all data
pct_cover <- data.frame(pct_cover=tapply(t_b_sums.df$t_part, list(t_b_sums.df$canopy), sum, na.rm=TRUE) / sum(sapply(d, function(i) nrow(i))) * 100)
write.csv(pct_cover, na='', file='pct_cover.csv')
print(pct_cover, digits=1)
write.csv(pct_cover, na='', file='pct_cover.csv')
print(pct_cover, digits=1)