R – Compute East-West or North-South Orientation of Polylines

coordinateslinestringrsf

I have a shapefile of road segments, which I want to segregate as East-West and North-South when viewing on Leaflet. The repex for the sf object is as follows:

enter image description here

dfr_sf <- structure(list(ID = c(1236, 761, 931, 1336, 762), geometry = structure(list(
    structure(c(-112.142904271646, -112.142904271255, -112.142904994711, 
    -112.142905247291, -112.142905736732, 33.5414077452841, 33.5414079101852, 
    33.5417621119912, 33.5418861392355, 33.5421262682246), .Dim = c(5L, 
    2L), class = c("XY", "LINESTRING", "sfg")), structure(c(-112.203748734908, 
    -112.203752166095, -112.203773879659, 33.4343700504272, 33.4342871568118, 
    33.433762597255), .Dim = 3:2, class = c("XY", "LINESTRING", 
    "sfg")), structure(c(-112.114883625126, -112.115167266906, 
    -112.115351914246, 33.6156119199386, 33.6158868152927, 33.6161422845864
    ), .Dim = 3:2, class = c("XY", "LINESTRING", "sfg")), structure(c(-111.964783456005, 
    -111.964090535608, -111.963936568584, -111.963468137801, 
    -111.962899435723, -111.962859452241, -111.962760305322, 
    -111.962671115879, -111.962595874294, 33.6187817091371, 33.6187867477904, 
    33.6187846321183, 33.6187781929798, 33.6187886003224, 33.6187893318816, 
    33.6187911449615, 33.6187927776802, 33.6187941550603), .Dim = c(9L, 
    2L), class = c("XY", "LINESTRING", "sfg")), structure(c(-112.121866995851, 
    -112.12203194047, -112.122267474906, -112.122538463309, -112.122772543993, 
    -112.122987457062, -112.123268503597, -112.123523950945, 
    -112.123670871747, -112.123825408105, -112.123972607689, 
    33.4804339199587, 33.4804340246712, 33.480434174243, 33.480434343499, 
    33.4804344897161, 33.4804346227382, 33.4804347963288, 33.4804349515442, 
    33.4804350401811, 33.4804351351541, 33.4804352257157), .Dim = c(11L, 
    2L), class = c("XY", "LINESTRING", "sfg"))), class = c("sfc_LINESTRING", 
"sfc"), precision = 0, bbox = structure(c(xmin = -112.203773879659, 
ymin = 33.433762597255, xmax = -111.962595874294, ymax = 33.6187941550603
), class = "bbox"), crs = structure(list(input = "+proj=longlat +datum=WGS84", 
    wkt = "GEOGCRS[\"unknown\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"longitude\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"latitude\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"), class = "crs"), n_empty = 0L)), row.names = c(827L, 
6013L, 387L, 1367L, 7384L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(ID = NA_integer_), .Label = c("constant", 
"aggregate", "identity"), class = "factor"))

I want to find each ID's or line's orientation in terms of either East-West or North-South and save the computation as EW or NS under a new variable in the sf object.

Updated based on comments

I approached this problem by identifying the range of latitudes and longitudes for each ID. If the range of longitudes is larger, it's a East-West orientation, and vice versa. Some IDs, however, had very small difference between the range of latitudes and longitudes, so their orientation could not be identified with certainty. Lines are roads from a grid-based system, so if some lines are radial, they could be labelled as EW or NS based on the angle (relative to north) between first and last points in the ID. If angle < 45 degree, it's NS, otherwise EW.

I'm hoping there's a better approach (perhaps in terms of both angle of polyline and ranges of lat/long) to figure out the orientation. Or is there a package-based function that identifies the orientation of a polyline based on its layout/orientation on a plot/map?

Best Answer

Some functions:

Compute the bearing using the geosphere package from the first to the last point of data defined as the first two columns in a matrix or data frame:

first_last_dir <- function(xy){
    geosphere::bearing(xy[1,1:2], xy[nrow(xy),1:2])
}

Get the coordinates of a spatial data frame of LINESTRING geometries, split by the L1 attribute (which breaks features), then pass to the above function to return the bearing:

orient <- function(lines){
    pts = data.frame(sf::st_coordinates(lines))
    pts = split(pts, pts$L1)
    bearing = sapply(pts, first_last_dir)
    bearing    
}

Then you can do:

dfr_sf$orient = orient(dfr_sf)

and make some plots with the bearing as title:

plots <- function(d){
    for(i in 1:nrow(d)){
        plot(d$geom[i])
        title(d$orient[i])
    }
}

enter image description here

Now code up a classification of N-S or E-W based on ranges of the numbers and you're done.

Note this uses geosphere::bearing which works give the initial bearing of the shortest distance path between the points based on the ellipsoid. Note the bearing from A to B might not be the same as from B to A because you are following something like a great circle, and the bearing changes as you go. Hence its possible (on this surface) that the bearing of A to B is E-W in your classification, but the bearing from B to A is N-S in the same classification. On a small scale this is unlikely to happen but its an unavoidable consequence of geometry on an ellipsoid.

Related Question