[R-sig-Geo] how to generate perpendicular transects along a line feature

[Barry Rowlingson]

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