r - 根据不同条件去除空间线段

Question

我正在尝试计算海岸线的暴露指数。我在沿海点周围每 5 度创建一条线。我已经擦除了与陆地相交的部分线条。但是，它会创建我不想要的线段。我需要：

(1) 如果它立即落入马达加斯加内陆，则排除整条线（红色）

(2) 选择线的第一段，例如删除在岛屿或任何陆地之后继续的线（绿色/蓝色）

（3）确保我与每条样线的点具有相同的ID（稍后我将针对许多点进行调整）

我无法选择特定的线段（即，如果线段与点具有相同的坐标）以及选择要删除的整条线。

见样带线和颜色参考

起始坐标

point
 <- class       : SpatialPointsDataFrame 
features    : 1 
extent      : 45.42639, 45.42639, -15.98098, -15.98098  (xmin, xmax, ymin, ymax)
crs         : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 5
names       :     layer,  path, Nearest_Sl, StdEr_SL, ID  
                                                                              

提取坐标和 ID

for (j in 1:length(point)){
  coords <- coordinates(point)
  ID <- point$ID
}

x <- cbind(ID, coords) 

从点计算线

library(sp)
library(geosphere)
library(spatstat)
library(maptools)

b=seq(0,355,5) # list bearings

# Calculate ending coordinate 
for(i in 1:length(b)){
  temp <- destPoint(p=coords,b=b[i],d=900000)# 900 km
  if(i==1){
    L <- cbind(x, temp)
  } else {
    L <- rbind(L,cbind(x, temp)) 
  }}

### Extracting beginning and end 
begin.coord <- data.frame(lon=c(L[,2]), lat=c(L[,3]))
end.coord <- data.frame(lon=c(L[,4]), lat=c(L[,5]))

### raw list to store Lines object
p <- psp(begin.coord[,1], begin.coord[,2], end.coord[,1], end.coord[,2],     owin(range(c(begin.coord[,1], end.coord[,1])), range(c(begin.coord[,2], end.coord[,2]))))

### Create spatial lines
p<- as(p, "SpatialLines")

### Remove line segments that overlap with world polygon
testclip <- raster::erase(p,world)

结果行的信息

testclip <-
class       : SpatialLines 
features    : 67 
extent      : 37.22043, 53.82955, -23.82036, -7.845263  (xmin, xmax, ymin, ymax)
crs         : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

### Example of 10th line with 6 segmented lines
str(testclip[10,])
Formal class 'SpatialLines' [package "sp"] with 3 slots
  ..@ lines      :List of 1
  .. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
  .. .. .. ..@ Lines:List of 6
  .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
  .. .. .. .. .. .. ..@ coords: num [1:2, 1:2] 45.4 48.8 -16 -12.6
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
  .. .. .. .. .. .. ..@ coords: num [1:2, 1:2] 48.8 48.8 -12.6 -12.6
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
  .. .. .. .. .. .. ..@ coords: num [1:2, 1:2] 48.9 49 -12.5 -12.4
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
  .. .. .. .. .. .. ..@ coords: num [1:2, 1:2] 49.1 49.2 -12.3 -12.2
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
  .. .. .. .. .. .. ..@ coords: num [1:2, 1:2] 49.2 49.2 -12.2 -12.1
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
  .. .. .. .. .. .. ..@ coords: num [1:2, 1:2] 49.3 51.2 -12.1 -10.2
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. ..@ ID   : chr "10"
  ..@ bbox       : num [1:2, 1:2] 45.4 -16 51.2 -10.2
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Testclip@lines[10]
[[1]]
An object of class "Lines"
Slot "Lines":
[[1]]
An object of class "Line"
Slot "coords":
            x         y
[1,] 45.42639 -15.98098
[2,] 48.82687 -12.56570


[[2]]
An object of class "Line"
Slot "coords":
            x         y
[1,] 48.83505 -12.55749
[2,] 48.83534 -12.55720


[[3]]
An object of class "Line"
Slot "coords":
            x         y
[1,] 48.89905 -12.49321
[2,] 48.95112 -12.44091


[[4]]
An object of class "Line"
Slot "coords":
            x         y
[1,] 49.12860 -12.26266
[2,] 49.15358 -12.23757


[[5]]
An object of class "Line"
Slot "coords":
            x         y
[1,] 49.23665 -12.15414
[2,] 49.24262 -12.14814


[[6]]
An object of class "Line"
Slot "coords":
            x         y
[1,] 49.33568 -12.05468
[2,] 51.22424 -10.15790



Slot "ID":
[1] "10"

 

标签： rspatialspspatstatgeosphere

Answer 1

首先，我会确保所有内容都投影到平面坐标（或使用地理纬度、经度坐标，以及最新版本的sf内置s2支持）。

如果您使用投影坐标，您可以spatstat大致使用如下：

1. 将海岸线及其周围的框视为线段的集合。
2. 从感兴趣的点画一条线，向一个方向延伸很远。
3. 检查方向是否在陆地上。
4. 如果不在陆地上，则找到延长线和海岸/箱体之间的所有交点。
5. 使用壁橱交点作为线的终点。
6. 重复 2.-5。对于每个感兴趣的角度。

library(spatstat)
# Some containing box:
box <- owin(c(330, 380), c(400, 450))
# Artificial landmass taken from spatstat dataset `chorley`:
landmass <- Window(chorley)
# Arbitrary point on coast:
pt1 <- midpoints.psp(edges(landmass)[1])
# Embed pt in big box:
Window(pt1) <- box
# Plot of setup:
plot(box)
plot(landmass, add = TRUE)
plot(pt1, add = TRUE, col = 2, cex = 1.5, pch = 20)

# Diameter of box for later usage:
D <- diameter(box)
# Numerical tolerance for later use:
eps <- 0.0001
# Box sides and coast as a single collection of line segments:
e <- superimpose(edges(box), edges(landmass))
# Angles to loop through:
angles <- seq(0, 350, by = 10)
# Object to hold ends for each angle:
ends <- ppp(window = box)
# Starting line from point extending far (`D`) towards East:
line0 <- as.psp(from = pt1, to = shift(pt1, vec = c(D,0)))
Window(line0) <- grow.rectangle(box, D)
# Test point just East of point:
test_pt0 <- shift(pt1, vec = c(eps, 0))
# Loop:
for(i in seq_along(angles)){
  # Rotate test point around coast point and check whether we are over land:
  test_pt <- rotate(test_pt0, angle = ang2rad(angles[i]), centre = pt1)
  overland <- inside.owin(test_pt, w = landmass)
  if(!overland){
    # Rotate starting line according to angle:
    line <- rotate(line0, angle = ang2rad(angles[i]), centre = pt1)
    # All crossing points between current line and edges of box+coast
    cross <- crossing.psp(line, e)
    # Index of closests crossing point (which is not the point it self):
    nn <- nncross(pt1, cross)
    if(nn$dist<eps){
      cross <- cross[-nn$which]
      nn <- nncross(pt1, cross)
    }
    # Add the point to the list of end points:
    ends <- superimpose(ends, cross[nn$which], W = box)
  }
}
# Lines from point to ends (only works if there are at least 2 end points:
rays1 <- as.psp(from = pt1, to = ends[1])
for(i in 2:npoints(ends)){
  rays1 <- superimpose(rays1, as.psp(from = pt1, to = ends[i]))
}
plot(e)
plot(pt1, add = TRUE, col = 2, cex = 1.5, pch = 20)
plot(rays1, add = TRUE, col = 4)

包装在一个函数中并应用于海岸上的一个新点：

rayfun <- function(pt){
  ends <- ppp(window = box)
  # Starting line from point extending far (`D`) towards East:
  line0 <- as.psp(from = pt, to = shift(pt, vec = c(D,0)))
  Window(line0) <- grow.rectangle(box, D)
  # Test point just East of point:
  test_pt0 <- shift(pt, vec = c(eps, 0))
  # Loop:
  for(i in seq_along(angles)){
    # Rotate test point around coast point and check whether we are over land:
    test_pt <- rotate(test_pt0, angle = ang2rad(angles[i]), centre = pt)
    overland <- inside.owin(test_pt, w = landmass)
    if(!overland){
      # Rotate starting line according to angle:
      line <- rotate(line0, angle = ang2rad(angles[i]), centre = pt)
      # All crossing points between current line and edges of box+coast
      cross <- crossing.psp(line, e)
      # Index of closests crossing point (which is not the point it self):
      nn <- nncross(pt, cross)
      if(nn$dist<eps){
        cross <- cross[-nn$which]
        nn <- nncross(pt, cross)
      }
      # Add the point to the list of end points:
      ends <- superimpose(ends, cross[nn$which], W = box)
    }
  }
  # Lines from point to ends (only works if there are at least 2 end points:
  rays <- as.psp(from = pt, to = ends[1])
  for(i in 2:npoints(ends)){
    rays <- superimpose(rays, as.psp(from = pt, to = ends[i]))
  }
  return(rays)
}
pt2 <- midpoints.psp(edges(landmass)[65])
rays2 <- rayfun(pt2)
plot(e)
plot(pt2, add = TRUE, col = 2, cex = 1.5, pch = 20)
plot(rays2, add = TRUE, col = 4)