Search code examples
rspatialr-spspatstatgeosphere

Remove segments of spatial lines based on different conditions


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

Starting coordinate

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  
                                                                              

Extracting coordinate and ID

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

x <- cbind(ID, coords) 

Calculate lines from point

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)

Information on resulting lines

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"

 

Solution

  • 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:

    1. Consider the coast line and the box around it as a collection of line segments.
    2. Make a line from the point of interest extending far i one direction.
    3. Check whether the direction is over land.
    4. If not over land find all intersection points between extended line and coast/box.
    5. Use the closet intersection point as end point of the line.
    6. Repeat 2.-5. for every angle of interest.
    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)