根据不同条件去除空间线段
Remove segments of spatial lines based on different conditions
我正在尝试计算海岸线的暴露指数。我在沿海点周围每隔 5 度创建了一条线。我已经擦掉了与陆地相交的部分线。但是，它创建了我不想要的线段。
我需要：
(1) 如果它立即落在马达加斯加内陆（红色），则排除整条线
(2) Select 该行的第一段
例如删除在岛屿或任何陆地之后继续的线 (green/blue)
(3) 确保我与每条样线的点具有相同的 ID
（后面我会调整很多点）
我无法选择特定的线段（即如果线段与点具有相同的坐标）以及选择要删除的整条线。
see transect lines and colour references
起始坐标
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"
首先我要确保一切都投影到平面坐标
（或对最新版本的所有内容使用地理纬度、经度坐标
sf 的 s2 支持 built-in)。
如果你使用投影坐标，你可以使用 spatstat 大致如下：
- 将海岸线及其周围的框视为线段的集合。
- 从兴趣点向一个方向延伸很远画一条线。
- 检查方向是否在陆地上。
- 如果不超过陆地，找到延长线和coast/box之间的所有交点。
- 使用壁橱交点作为线的终点。
- 重复 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)
我正在尝试计算海岸线的暴露指数。我在沿海点周围每隔 5 度创建了一条线。我已经擦掉了与陆地相交的部分线。但是，它创建了我不想要的线段。 我需要：
(1) 如果它立即落在马达加斯加内陆（红色），则排除整条线
(2) Select 该行的第一段 例如删除在岛屿或任何陆地之后继续的线 (green/blue)
(3) 确保我与每条样线的点具有相同的 ID （后面我会调整很多点）
我无法选择特定的线段（即如果线段与点具有相同的坐标）以及选择要删除的整条线。
see transect lines and colour references
起始坐标
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"
首先我要确保一切都投影到平面坐标
（或对最新版本的所有内容使用地理纬度、经度坐标
sf 的 s2 支持 built-in)。
如果你使用投影坐标，你可以使用 spatstat 大致如下：
- 将海岸线及其周围的框视为线段的集合。
- 从兴趣点向一个方向延伸很远画一条线。
- 检查方向是否在陆地上。
- 如果不超过陆地，找到延长线和coast/box之间的所有交点。
- 使用壁橱交点作为线的终点。
- 重复 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)