Determine transects perpendicular to a (coast)line in R

Question

I'd like to automatically derive transects, perpendicular to the coastline. I need to be able to control their length and spacing and their oriëntation needs to be on the "correct" side of the line. I came up with a way to do that, but especially selecting the "correct" (it needs to point to the ocean) can be done better. General approach:

1. For each line segment in a SpatialLineDataFrame define transect locations
2. define transect: in both directions perpendicular to coastline: create points that determine the transect
3. Create a polygon based on the coastline, add extra points to grow the polygon in a direction that is known and use that to clip the points that are inside (considered as land, and therefore not of interest)
4. Create transect based on remaining point

Especially part 3 is of interest. I'd like a more robust method to determine the correct orientation of the transect. This is what i'm using now:

    library(rgdal)
library(raster)
library(sf)
library(ggplot2)
library(rgeos)      # create lines and spatial objects

# create testing lines
l1 <- cbind(c(1, 2, 3), c(3, 2, 2))
l2 <- cbind(c(1, 2, 3), c(1, 1.5, 1))

Sl1 <- Line(l1)
Sl2 <- Line(l2)

S1 <- Lines(list(Sl1), ID = "a")
S2 <- Lines(list(Sl2), ID = "b")

line <- SpatialLines(list(S1, S2))
plot(line)

# for testing:
sep <- 0.1
start <- 0


AllTransects <- vector('list', 100000) # DB that should contain all transects

for (i in 1: length(line)){
    # i <- 2

    ###### Define transect locations
    # Define geometry subset
    subset_geometry <- data.frame(geom(line[i,]))[, c('x', 'y')]

    # plot(SpatialPoints(data.frame(x = subset_geometry[,'x'], y = subset_geometry[,'y'])), axes = T, add = T)

    dx <- c(0, diff(subset_geometry[,'x'])) # Calculate difference at each cell comapred to next cell
    dy <- c(0, diff(subset_geometry[,'y']))

    dseg <- sqrt(dx^2+dy^2)                 # get rid of negatives and transfer to uniform distance per segment (pythagoras)
    dtotal <- cumsum(dseg)                  # cumulative sum total distance of segments

    linelength = sum(dseg)                  # total linelength
    pos = seq(start,linelength, by=sep)     # Array with postions numbers in meters
    whichseg = unlist(lapply(pos, function(x){sum(dtotal<=x)})) # Segments corresponding to distance

    pos=data.frame(pos=pos,                            # keep only 
                   whichseg=whichseg,                  # Position in meters on line
                   x0=subset_geometry[whichseg,1],     # x-coordinate on line
                   y0=subset_geometry[whichseg,2],     # y-coordinate on line
                   dseg = dseg[whichseg+1],            # segment length selected (sum of all dseg in that segment)
                   dtotal = dtotal[whichseg],          # Accumulated length
                   x1=subset_geometry[whichseg+1,1],   # Get X coordinate on line for next point
                   y1=subset_geometry[whichseg+1,2]    # Get Y coordinate on line for next point
    )

    pos$further =  pos$pos - pos$dtotal       # which is the next position (in meters)
    pos$f = pos$further/pos$dseg              # fraction next segment of its distance
    pos$x = pos$x0 + pos$f * (pos$x1-pos$x0)  # X Position of point on line which is x meters away from x0
    pos$y = pos$y0 + pos$f * (pos$y1-pos$y0)  # Y Position of point on line which is x meters away from y0

    pos$theta = atan2(pos$y0-pos$y1,pos$x0-pos$x1)  # Angle between points on the line in radians
    pos$object = i

    ###### Define transects
    tlen <- 0.5
    pos$thetaT = pos$theta+pi/2         # Get the angle
    dx_poi <- tlen*cos(pos$thetaT) # coordinates of point of interest as defined by position length (sep)
    dy_poi <- tlen*sin(pos$thetaT) 

    # tabel met alleen de POI informatie
    # transect is defined by x0,y0 and x1,y1 with x,y the coordinate on the line
    output <-     data.frame(pos = pos$pos,
                             x0 = pos$x + dx_poi,       # X coordinate away from line
                             y0 = pos$y + dy_poi,       # Y coordinate away from line
                             x1 = pos$x - dx_poi,       # X coordinate away from line
                             y1 = pos$y - dy_poi,       # X coordinate away from line
                             theta = pos$thetaT,    # angle
                             x = pos$x,             # Line coordinate X
                             y = pos$y,             # Line coordinate Y
                             object = pos$object,
                             nextx = pos$x1,
                             nexty = pos$y1) 

    # create polygon from object to select correct segment of the transect (coastal side only) 
    points_for_polygon <- rbind(output[,c('x', 'y','nextx', 'nexty')])# select points
    pol_for_intersect <- SpatialPolygons( list( Polygons(list(Polygon(points_for_polygon[,1:2])),1)))
    # plot(pol_for_intersect, axes = T, add = T)

    # Find a way to increase the polygon - should depend on the shape&direction of the polygon
    # for the purpose of cropping the transects
    firstForPlot <- data.frame(x = points_for_polygon$x[1], y = points_for_polygon$y[1])
    lastForPlot <- data.frame(x = points_for_polygon$x[length(points_for_polygon$x)],
                              y = points_for_polygon$y[length(points_for_polygon$y)])

    plot_first <- SpatialPoints(firstForPlot)
    plot_last <- SpatialPoints(lastForPlot)
    # plot(plot_first, add = T, col = 'red')
    # plot(plot_last, add = T, col = 'blue')

    ## Corners of shape dependent bounding box
    ## absolute values should be depended on the shape beginning and end point relative to each other??
    LX <- min(subset_geometry$x)
    UX <- max(subset_geometry$x)
    LY <- min(subset_geometry$y)
    UY <- max(subset_geometry$y)
    # polygon(x = c(LX, UX, UX, LX), y = c(LY, LY, UY, UY), lty = 2)
    # polygon(x = c(LX, UX, LX), y = c(LY, LY, UY), lty = 2)

    # if corners are changed to much the plot$near becomes a problem: the new points are to far away
    # Different points are selected
    LL_corner <- data.frame(x = LX-0.5, y = LY - 1)
    LR_corner <- data.frame(x = UX + 0.5 , y = LY - 1)
    UR_corner <- data.frame(x = LX, y = UY)
    corners <- rbind(LL_corner, LR_corner)
    bbox_add <- SpatialPoints(rbind(LL_corner, LR_corner))
    # plot(bbox_add ,col = 'green', axes = T, add = T)

    # Select nearest point for drawing order to avoid weird shapes
    firstForPlot$near <-apply(gDistance(bbox_add,plot_last, byid = T), 1, which.min)
    lastForPlot$near <- apply(gDistance(bbox_add,plot_first, byid = T), 1, which.min)

    # increase polygon with corresponding points
    points_for_polygon_incr <- rbind(points_for_polygon[1:2], corners[firstForPlot$near,], corners[lastForPlot$near,])
    pol_for_intersect_incr <- SpatialPolygons( list( Polygons(list(Polygon(points_for_polygon_incr)),1)))
    plot(pol_for_intersect_incr, col = 'blue', axes = T)

    # Coordinates of points first side
    coordsx1y1 <- data.frame(x = output$x1, y = output$y1)
    plotx1y1 <- SpatialPoints(coordsx1y1)
    plot(plotx1y1, add = T)

    coordsx0y0 <- data.frame(x = output$x0, y = output$y0)
    plotx0y0 <- SpatialPoints(coordsx0y0)
    plot(plotx0y0, add = T, col = 'red')

    # Intersect
    output[, "x1y1"] <- over(plotx1y1, pol_for_intersect_incr)
    output[, "x0y0"] <- over(plotx0y0, pol_for_intersect_incr)   
    x1y1NA <- sum(is.na(output$x1y1)) # Count Na  
    x0y0NA <- sum(is.na(output$x1y1)) # Count NA

    # inefficient way of selecting the correct end point
    # e.g. either left or right, depending on intersect
    indexx0y0 <- with(output, !is.na(output$x0y0))
    output[indexx0y0, 'endx'] <- output[indexx0y0, 'x1']
    output[indexx0y0, 'endy'] <- output[indexx0y0, 'y1']

    index <- with(output, is.na(output$x0y0))
    output[index, 'endx'] <- output[index, 'x0']
    output[index, 'endy'] <- output[index, 'y0']

    AllTransects = rbind(AllTransects, output)
}



# Create the transects
lines <- vector('list', nrow(AllTransects))
for(n in 1: nrow(AllTransects)){
  # n = 30

  begin_coords <- data.frame(lon = AllTransects$x, lat = AllTransects$y)       # Coordinates on the original line
  end_coords <- data.frame(lon = AllTransects$endx, lat = AllTransects$endy)   # coordinates as determined by the over: remove implement in row below by selecting correct column from output

  col_names <- list('lon', 'lat')
  row_names <- list('begin', 'end')
  # dimnames < list(row_names, col_names)

  x <- as.matrix(rbind(begin_coords[n,], end_coords[n,]))

  dimnames(x) <- list(row_names, col_names)
  lines[[n]] <- Lines(list(Line(x)), ID = as.character(n))
}
lines_sf <- SpatialLines(lines)
# plot(lines_sf)
df <- SpatialLinesDataFrame(lines_sf, data.frame(AllTransects))

plot(df, axes = T)

As long as i'm able to correctly define the bounding box and grow the polygon correctly this works. But I'd like to try this on multiple coastlines and parts of coastlines, each with its own orientation. In the example below the growing of the polygon is made for the bottom coastline segment, as a result the top one has transects in the wrong direction.

Anybody has an idea in what directio to look? I was considering to perhaps use external data but when possible i'd like to avoid that.

Answer 1

I used your code for my question (measure line inside a polygon) but maybe this works for you:

1. Took a spatial polygon or line
2. Extract the coordinates of the element
3. Make a combination of coordinates to create straight lines, from with you can derivate perpendicular lines (e.g. ((x1,x3)(y1, y3)) or ((x2,x4)(y2, y4)) )
4. Iterate along with all the pairs of coordinates
5. Apply the code you did, especially the result of the 'output' table.

I did this for a polygon, so I could generate perpendicular lines based on the straight line I create taking an arbitrary (1, 3) set of coordinates.

#Define a polygon
pol <- rip[1, 1] # I took the first polygon from my Shapefile
polcoords <- pol@polygons[[1]]@Polygons[[1]]@coords

# define how to create your coords pairing. My case: 1st with 3rd, 2nd with 4th, ...
pairs <- data.frame(a = 1:( nrow(polcoords) - 1), 
                    b = c(2:(nrow(polcoords)-1)+1, 1) )

# Empty list to store the lines
lnDfls <- list()

for (j in 1:nrow(pairs)){ # j = 1
  
  # Select the pairs
  pp <- polcoords[c(pairs$a[j], pairs$b[j]), ]
  
  #Extract mean coord, from where the perp. line will start
  midpt <- apply(pp, 2, mean)
  
  # points(pp, col = 3, pch = 20 )
  # points(midpt[1], midpt[2], col = 4, pch = 20)
  
  x <- midpt[1]
  y <- midpt[2]

  theta = atan2(y = pp[2, 2] - pp[1, 2], pp[2, 1] - pp[1, 1])  # Angle between points on the line in radians
  # pos$theta = atan2(y = pos$y0-pos$y1 , pos$x0-pos$x1)  # Angle between points on the line in radians
  
  ###### Define transects
  tlen <- 1000 # distance in m
  thetaT = theta+pi/2         # Get the angle
  dx_poi <- tlen*cos(thetaT) # coordinates of point of interest as defined by position length (sep)
  dy_poi <- tlen*sin(thetaT) 
  
  # tabel met alleen de POI informatie
  # transect is defined by x0,y0 and x1,y1 with x,y the coordinate on the line
  output2 <-     data.frame(#pos = pos,
                           x0 = x + dx_poi,       # X coordinate away from line
                           y0 = y + dy_poi,       # Y coordinate away from line
                           x1 = x - dx_poi,       # X coordinate away from line
                           y1 = y - dy_poi       # X coordinate away from line
                           #theta = thetaT,    # angle
                           #x = x,             # Line coordinate X
                           #y = y             # Line coordinate Y
                           ) 
  
  # points(output2$x1, output2$y1, col = 2)
  #segments(x, y, output2$x1[1], output2$y1[1], col = 2)
  
  mat <- as.matrix(cbind( c( x, output2$x1[1] ) , c( y, output2$y1[1] ) ))
  
  LL <- Lines(list(Line( mat )), ID = as.character(j))
  # plot(SpatialLinesDataFrame(LL, data.frame (a = 1)), add = TRUE, col = 2)
  # plot(SpatialLines(list(LL)), add = TRUE, col = 2)
  
  
  #lnList[[j]] <- LL
  lnDfls[[j]] <- SpatialLinesDataFrame( SpatialLines(LinesList = list(LL)) ,
                                        match.ID = FALSE,
                                      data.frame(id = as.character(j ) ) )
   
  # line = st_sfc(st_linestring(mat))
  # st_length(line)
   
  # ln <- (SpatialLines(LinesList = list(LL)))
  # lndf <- SpatialLinesDataFrame(  lndf , data.frame(id = j ))
  # sf::st_length(ln)
  # # plot(lines_sf)
}

compDf <- do.call(what = sp::rbind.SpatialLines, args = lnDfls)


plot(pol)
plot(compDf, add = TRUE, col = 2)
plot(inDfLn, add = TRUE, col  = 3)