问题描述
我正在尝试计算海岸线的暴露指数。我围绕沿海点每5度创建了一条线。我已删除了与土地相交的部分线。但是，它会创建我不想要的线段。 我需要：
（1）如果整条线立即落入马达加斯加内陆（红色），则将其排除在外
（2）选择该行的第一段 例如删除在岛屿或任何土地（绿色/蓝色）之后继续延伸的线段
（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
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]),5]))
### raw list to store Lines object
p <- psp(begin.coord[,1],begin.coord[,2],end.coord[,owin(range(c(begin.coord[,1])),range(c(begin.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,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(*,1:2] 48.9 49 -12.5 -12.4
.. .. .. .. .. .. .. ..- attr(*,1:2] 49.1 49.2 -12.3 -12.2
.. .. .. .. .. .. .. ..- attr(*,1:2] 49.2 49.2 -12.2 -12.1
.. .. .. .. .. .. .. ..- attr(*,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"
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中的一个。
如果使用投影坐标，则可以大致如下使用spatstat：
- 将海岸线及其周围的框视为线段的集合。
- 从兴趣点开始沿一条方向延伸一条线。
- 检查方向是否在陆地上。
- 如果不在陆地上，请找到延长线与海岸/盒子之间的所有交点。
- 使用壁橱相交点作为直线的终点。
- 重复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,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,pch = 20)
plot(rays1,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,0)))
Window(line0) <- grow.rectangle(box,D)
# Test point just East of point:
test_pt0 <- shift(pt,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,centre = pt)
overland <- inside.owin(test_pt,w = landmass)
if(!overland){
# Rotate starting line according to angle:
line <- rotate(line0,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,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,pch = 20)
plot(rays2,col = 4)
