Search code examples
rfiltergpsregion

R filter GPS dataset by vector with margins, or a bounding box at an angle


I have a large dataset which includes longitude and latitude values which I'm analysing in R. e.g.

structure(list(S = structure(c(1L, 1L, 1L, 1L, 1L, 1L), .Label = "L", class = "factor"), 
No = 7:12, X.. = c(175L, 200L, 225L, 250L, 275L, 300L), X. = c(200L, 
225L, 250L, 275L, 300L, 325L), Time = c(112836.789, 112837.291, 
112837.793, 112838.295, 112838.797, 112839.299), Latidude = c(51.28769291, 
51.28768926, 51.28768518, 51.28768096, 51.28767662, 51.28767237
), Longitude = c(0.868181943, 0.868182069, 0.868181418, 0.868179761, 
0.868177937, 0.868176734), Altitude = c(61.177, 61.145, 61.104, 
61.06, 61.05, 61.03), Bearing = c(201.3103, 201.3103, 201.3103, 
201.3103, 201.3103, 201.3103), Tree = structure(c(1L, 1L, 
1L, 1L, 1L, 1L), .Label = "Tree", class = "factor"), XSArea = c(1.4204362, 
0.4883868, 0.5491678, 0.6008296, 0.5633923, 0.5917539), AveWidth = c(1.1840309, 
0.5247761, 0.7371026, 1.2021244, 0.7507043, 0.3796109), AreaDensity = c(0.1989031, 
3.6698774, 3.7629885, 2.8899686, 2.5806617, 1.5068812)), row.names = c(NA,6L), class = "data.frame") 

I'm trying to filter this dataset to a small region, which is a rectangle at an angle. I have a start point and an end point, plus a bearing and rectangle width. I don't know how to create the bounding box at an angle or then filter the large dataset by this bounding box. An example of the bounding box is:

bb<-cbind(c(0.86777,51.28743),c(0.86744,51.28690))
require(ggmaps)
bounding_box<-make_bbox(bb[1,],bb[2,])

But that just makes a rectangle running North to South, and i need it at an angle

require(geosphere)
bearingRhumb(bb[,1], bb[,2])

200.5686 So the bounding box rectangle should be at that angle

I've tried a different approach by creating a vector of longitudes and latitudes, which represent a straight line.

vect_line<-gcIntermediate(bb[,1], bb[,2], n = 1000, addStartEnd = TRUE)

This straight line is calculated from two GPS points using 'gcIntermediate' from the geosphere package. I planned to then filter the large dataset by this straight line, but i need to add a margin to the straight line. So that instead of the filter being limited to matching exactly each value in the vector, it filters by the value in the vector plus and minus a bit. The plus or minus a bit is a distance (e.g. 1 m) around each of the GPS coordinates within my lon and lat vector. Perhaps something like a fuzzy filter? Or maybe expand the straight line vector somehow?

I'm open to any ideas really, so if anyone has a better way to do what I'm trying to do that is also fine. e.g. create a bounding box, which is a rectangle at an angle (the bearing of my straight line), and filter my large dataset using the bounding box.

Let me know if i need to add more details to my question.


Solution

  • From your question, I assume your problem lies with creating the box you want to filter your large dataset on...

    So this answer focusses of creating this bbox... It does not do the actual filtering.. please let me know if you cannot get this done yourself.

    here we go..

    library( geosphere )
    library( sf )
    #boundary box sample coordinates
    bb1 <- c( 0.86777, 51.28743 )
    bb2 <- c( 0.86744, 51.28690 )
    #set total width (in meters) of the boundary-box
    bb.width = 50
    #calculate course from 1st >> 2nd bb point
    bb.bearing <- ( geosphere::bearing( bb1, bb2 ) + 360 ) %% 360
    #perpendicular angles
    bb.bearing.perp1 <- (bb.bearing + 360 + 90) %% 360
    bb.bearing.perp2 <- (bb.bearing + 360 - 90) %% 360
    #calculate the four corner-coordinates of the boundary-box
    # leave form the start- and end point of a line, and travel
    # half the desired width of the box in a perpendicular direction
    bb.c1 <- geosphere::destPoint( p = bb1, bb.bearing.perp1, 0.5 * bb.width )
    bb.c2 <- geosphere::destPoint( p = bb1, bb.bearing.perp2, 0.5 * bb.width )
    bb.c3 <- geosphere::destPoint( p = bb2, bb.bearing.perp1, 0.5 * bb.width )
    bb.c4 <- geosphere::destPoint( p = bb2, bb.bearing.perp2, 0.5 * bb.width )
    #make spatial object out of the four points, ans make it a real 'box'
    bb.sf <- rbind( bb.c1, bb.c2, bb.c3, bb.c4 ) %>%
      as.data.frame() %>%
      sf::st_as_sf( coords = c("lon", "lat"), crs = 4326 ) %>%
      sf::st_union() %>% 
      sf::st_convex_hull()
    

    now, perform a visual check of what we have done

    #view original bb.points bb1 and bb2, for viewing on the map
    bb.points.sf <- rbind( bb1, bb2 ) %>% 
      as.data.frame() %>% 
      st_as_sf( coords = c("V1", "V2"), crs = 4326 )
    
    mapview::mapview( list( bb.points.sf, bb.sf ) )
    

    enter image description here

    looks pretty solid to me ;-)