Определить разрезы перпендикулярно (береговой) линии в R

Question

Я бы хотел автоматически получать трансекты, перпендикулярные береговой линии. Мне нужно иметь возможность контролировать их длину и расстояние, а их ориентация должна быть на «правильной» стороне линии. Я придумал способ сделать это, но особенно выбор «правильный» (он должен указывать на океан) может быть сделан лучше. Общий подход:

1. Для каждого сегмента линии в SpatialLineDataFrame определите трансект Места
2. определить трансект: в обоих направлениях, перпендикулярных береговой линии: создать точки, определяющие трансект
3. Создание многоугольника на основе береговой линии, добавление дополнительных точек для роста многоугольника в известном направлении и использование его для обрезки точек, которые находятся внутри (считаются землей и поэтому не представляют интереса)
4. Создать трансект на основе оставшейся точки

Особенно интересна часть 3. Я бы хотел более надежный метод определения правильной ориентации разреза. Это то, что я сейчас использую:

    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)

Пока я могу правильно определить ограничивающую рамку и правильно вырастить многоугольник, это работает. Но я хотел бы попробовать это на нескольких береговых линиях и участках береговых линий, каждая со своей ориентацией. В приведенном ниже примере выращивание многоугольника выполняется для нижнего сегмента береговой линии, в результате чего верхний имеет пересечения в неправильном направлении.

У кого-нибудь есть идея, в каком направлении смотреть? Я подумывал, возможно, использовать внешние данные, но, когда это возможно, я бы хотел этого избежать.