[R-sig-Geo] how to generate perpendicular transects along a line feature
Barry Rowlingson
b.rowlingson at lancaster.ac.uk
Wed Jun 5 14:36:28 CEST 2013
Nothing else to do this lunchtime...
Two functions - one computes evenly spaced points along a set of xy
coordinates, the other works out the end point of transects centred on
those points:
evenspace <- function(xy, sep, start=0, size){
dx <- c(0,diff(xy[,1]))
dy <- c(0,diff(xy[,2]))
dseg <- sqrt(dx^2+dy^2)
dtotal <- cumsum(dseg)
linelength = sum(dseg)
pos = seq(start,linelength, by=sep)
whichseg = unlist(lapply(pos, function(x){sum(dtotal<=x)}))
pos=data.frame(pos=pos,whichseg=whichseg,
x0=xy[whichseg,1],
y0=xy[whichseg,2],
dseg = dseg[whichseg+1],
dtotal = dtotal[whichseg],
x1=xy[whichseg+1,1],
y1=xy[whichseg+1,2]
)
pos$further = pos$pos - pos$dtotal
pos$f = pos$further/pos$dseg
pos$x = pos$x0 + pos$f * (pos$x1-pos$x0)
pos$y = pos$y0 + pos$f * (pos$y1-pos$y0)
pos$theta = atan2(pos$y0-pos$y1,pos$x0-pos$x1)
return(pos[,c("x","y","x0","y0","x1","y1","theta")])
}
transect <- function( tpts, tlen){
tpts$thetaT = tpts$theta+pi/2
dx = tlen*cos(tpts$thetaT)
dy = tlen*sin(tpts$thetaT)
return(
data.frame(x0 = tpts$x + dx,
y0 = tpts$y + dy,
x1 = tpts$x - dx,
y1 = tpts$y -dy)
)
}
Apologies for complete lack of documentation and comments! If I had
two lunchtimes this would be all written up, documented, and packaged
onto GitHub....
Sample usage:
1. make a meandering stream:
> stream = data.frame(x=seq(0,20,len=50))
> stream$y=sin(stream$x/3)
2. plot it - note asp=1 is essential:
> plot(stream,type="l", asp=1)
3. compute evenly spaced (1.5 units) transect centres:
> tspts = evenspace(stream,1.5)
> points(tspts,pch="X")
4. compute perpendicular transects of length 2:
> tslines = transect(tspts,2)
> segments(tslines$x0,tslines$y0,tslines$x1,tslines$y1)
the transect function takes the output of the evenspace function and
returns a data frame with the end points suitable for feeding into
segments(). You can convert all these things into sp-class
SpatialLines objects and save them as Shapefiles...
Nice example, I might write this up as a blog post on that RStudio
blog thing I've forgotten the name of, or use it for some teaching I'm
doing.
The next interesting step would be to use rgeos:gBuffer to generate
buffer zones round the transects....
Barry
More information about the R-sig-Geo
mailing list