I am trying to calculate an exposure index of the coastline. I have created lines every 5 degrees around a coastal point. I have erased the parts of the lines that intersect over land. However, it creates line segments that I do not want. I need to:
(1) Exclude the whole line if it immediately falls inland of Madagascar (red)
(2) Select the first segment of the line e.g. Remove lines that continue after an island or any land (green/blue)
(3) Make sure I have the same ID as the point for each transect line (later I will adjust this for many points)
I am having trouble selecting specific segments of lines (i.e. if the segment of the line has the same coordinates as a point) as well as selecting the whole line to remove.
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
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"
First of all I would make sure everything is projected to planar coordinates
(or use geographic lat,lon coordinates for everything with a recent version
of sf
with s2
support built-in).
If you go with projected coordinates you can use spatstat
roughly as follows:
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)
Wrapped in a function and applied to a new point on the coast:
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)