Find nearest line segment to each point in R

Question

I've seen similar questions, but not with R. I have a sf object with approx 500000 thousand points. I also have an sf object of lines. Each line has a unique segment number called 'segnum'. For each pixel, I want to know which line (identified by segment number) is the closest. The point and line files are currently in geographic coordinates (crs = 4326), but I can reproject them if need be.

Based on this question (How to calculate distance from POINT to LINESTRING in R using `sf` library and get all POINT features from LINESTRING where distances were calculated?) I was thinking of using dist2line from geosphere.

For example: dist = geosphere::dist2Line(p=st_coordinates(pixels), line = st_coordinates(lines)[,1:2])

But that just gives you the distance from the point to the line. It doesn't tell you which line segment is closest. Any ideas of how to do this?

Here is a small sample of my data:

pixels = structure(list(id = c("1", "1", "1", "1", "1"), system.index = 0:4, 
    b2_allPix = c(7009L, 3932L, 7329L, 5734L, 6525L), b2_cloudfilt = c(0L, 
    3932L, 0L, 0L, 0L), system.time_start = c(951350400000, 951350400000, 
    951350400000, 951350400000, 951350400000), .geo = c(NA, NA, 
    NA, NA, NA), geometry = structure(list(`1` = structure(c(-137.849125851618, 
    66.7968749939405), class = c("XY", "POINT", "sfg")), `2` = structure(c(-137.838550341857, 
    66.7968749939405), class = c("XY", "POINT", "sfg")), `3` = structure(c(-137.879732491708, 
    66.7947916606074), class = c("XY", "POINT", "sfg")), `4` = structure(c(-137.862752400773, 
    66.7927083272742), class = c("XY", "POINT", "sfg")), `5` = structure(c(-137.881661018612, 
    66.7885416606079), class = c("XY", "POINT", "sfg"))), .Names = c("1", 
    "2", "3", "4", "5"), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(-137.881661018612, 
    66.7885416606079, -137.838550341857, 66.7968749939405), .Names = c("xmin", 
    "ymin", "xmax", "ymax"), class = "bbox"), crs = structure(list(
        epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), .Names = c("epsg", 
    "proj4string"), class = "crs"), n_empty = 0L)), .Names = c("id", 
"system.index", "b2_allPix", "b2_cloudfilt", "system.time_start", 
".geo", "geometry"), row.names = c(NA, 5L), class = c("sf", "data.frame"
), sf_column = "geometry", agr = structure(c(1L, 1L, 1L, 1L, 
1L, 1L), .Names = c("id", "system.index", "b2_allPix", "b2_cloudfilt", 
"system.time_start", ".geo"), .Label = c("constant", "aggregate", 
"identity"), class = "factor")) 

lines = structure(list(id = structure(c(1L, 112L, 223L, 334L, 445L), .Label = c("1", 
"10", "100", "101", "102", "103", "104", "105", "106", "107", 
"108", "109", "11", "110", "111", "112", "113", "114", "115", 
"116", "117", "118", "119", "12", "120", "121", "122", "123", 
"124", "125", "126", "127", "128", "129", "13", "130", "131", 
"132", "133", "134", "135", "136", "137", "138", "139", "14", 
"140", "141", "142", "143", "144", "145", "146", "147", "148", 
"149", "15", "150", "151", "152", "153", "154", "155", "156", 
"157", "158", "159", "16", "160", "161", "162", "163", "164", 
"165", "166", "167", "168", "169", "17", "170", "171", "172", 
"173", "174", "175", "176", "177", "178", "179", "18", "180", 
"181", "182", "183", "184", "185", "186", "187", "188", "189", 
"19", "190", "191", "192", "193", "194", "195", "196", "197", 
"198", "199", "2", "20", "200", "201", "202", "203", "204", "205", 
"206", "207", "208", "209", "21", "210", "211", "212", "213", 
"214", "215", "216", "217", "218", "219", "22", "220", "221", 
"222", "223", "224", "225", "226", "227", "228", "229", "23", 
"230", "231", "232", "233", "234", "235", "236", "237", "238", 
"239", "24", "240", "241", "242", "243", "244", "245", "246", 
"247", "248", "249", "25", "250", "251", "252", "253", "254", 
"255", "256", "257", "258", "259", "26", "260", "261", "262", 
"263", "264", "265", "266", "267", "268", "269", "27", "270", 
"271", "272", "273", "274", "275", "276", "277", "278", "279", 
"28", "280", "281", "282", "283", "284", "285", "286", "287", 
"288", "289", "29", "290", "291", "292", "293", "294", "295", 
"296", "297", "298", "299", "3", "30", "300", "301", "302", "303", 
"304", "305", "306", "307", "308", "309", "31", "310", "311", 
"312", "313", "314", "315", "316", "317", "318", "319", "32", 
"320", "321", "322", "323", "324", "325", "326", "327", "328", 
"329", "33", "330", "331", "332", "333", "334", "335", "336", 
"337", "338", "339", "34", "340", "341", "342", "343", "344", 
"345", "346", "347", "348", "349", "35", "350", "351", "352", 
"353", "354", "355", "356", "357", "358", "359", "36", "360", 
"361", "362", "363", "364", "365", "366", "367", "368", "369", 
"37", "370", "371", "372", "373", "374", "375", "376", "377", 
"378", "379", "38", "380", "381", "382", "383", "384", "385", 
"386", "387", "388", "389", "39", "390", "391", "392", "393", 
"394", "395", "396", "397", "398", "399", "4", "40", "400", "401", 
"402", "403", "404", "405", "406", "407", "408", "409", "41", 
"410", "411", "412", "413", "414", "415", "416", "417", "418", 
"419", "42", "420", "421", "422", "423", "424", "425", "426", 
"427", "428", "429", "43", "430", "431", "432", "433", "434", 
"435", "436", "437", "438", "439", "44", "440", "441", "442", 
"443", "444", "445", "446", "447", "448", "449", "45", "450", 
"451", "452", "453", "454", "455", "456", "457", "458", "459", 
"46", "460", "461", "462", "463", "464", "465", "466", "467", 
"468", "469", "47", "470", "471", "472", "473", "474", "475", 
"476", "477", "478", "479", "48", "480", "481", "482", "483", 
"484", "485", "486", "487", "488", "489", "49", "490", "491", 
"492", "493", "494", "495", "496", "497", "498", "499", "5", 
"50", "500", "501", "502", "503", "504", "505", "506", "507", 
"508", "509", "51", "510", "52", "53", "54", "55", "56", "57", 
"58", "59", "6", "60", "61", "62", "63", "64", "65", "66", "67", 
"68", "69", "7", "70", "71", "72", "73", "74", "75", "76", "77", 
"78", "79", "8", "80", "81", "82", "83", "84", "85", "86", "87", 
"88", "89", "9", "90", "91", "92", "93", "94", "95", "96", "97", 
"98", "99"), class = "factor"), segnum = 1:5, geometry = structure(list(
    structure(c(-163.939480000102, -163.932950000242, -163.930539999721, 
    -163.93100000025, -163.933320000218, -163.935640000186, -163.941510000186, 
    -163.947380000187, -163.950175000026, -163.952969999866, 
    -163.952909999797, -163.948820000043, -163.953614999932, 
    -163.958409999822, -163.958460000329, -163.955969999716, 
    -163.948809999582, -163.950620000313, -163.94980000027, -163.945235000195, 
    -163.94067000012, -163.93931771722, 62.4387199999281, 62.4343600003138, 
    62.4286899996438, 62.4238500003773, 62.4198250000238, 62.4157999996703, 
    62.4099199996583, 62.4040399996463, 62.3988099999316, 62.3935800002169, 
    62.387119999988, 62.3795499998327, 62.3751399997113, 62.3707299995899, 
    62.3691199999881, 62.3658599998392, 62.3633900001485, 62.3612500003873, 
    62.3598899997237, 62.3578349998354, 62.3557799999471, 62.3541680789363
    ), .Dim = c(22L, 2L), class = c("XY", "LINESTRING", "sfg"
    )), structure(c(-163.93931771722, -163.938169999946, -163.93740000041, 
    -163.938740000151, -163.944244999957, -163.949749999763, 
    -163.947819999793, -163.939500000125, -163.945110000126, 
    -163.950720000128, -163.956330000129, -163.959765000255, 
    -163.96320000038, -163.967020000273, -163.970840000166, -163.974660000059, 
    -163.978479999952, -163.979979999877, -163.981479999802, 
    -163.985049999858, -163.988619999913, -163.988893564025, 
    62.3541680789363, 62.3528000001199, 62.3498299998546, 62.3447300002894, 
    62.3405950002591, 62.3364600002289, 62.3337599998242, 62.3315399999711, 
    62.3283500000525, 62.325160000134, 62.3219700002154, 62.3189000002849, 
    62.3158300003543, 62.3098650002446, 62.303900000135, 62.2979350000253, 
    62.2919699999156, 62.2861900000185, 62.2804100001214, 62.2773400001909, 
    62.2742700002603, 62.2739300405841), .Dim = c(22L, 2L), class = c("XY", 
    "LINESTRING", "sfg")), structure(c(-163.988893564025, -163.992389999749, 
    -163.996159999584, -163.996669999721, -163.994829999855, 
    -163.992989999988, -163.994169999995, -163.995350000002, 
    -163.999320000068, -164.003290000133, -164.005879999962, 
    -164.005639999686, -164.001679999632, -163.995054999888, 
    -163.988430000144, -163.980987500023, -163.973544999902, 
    -163.966102499781, -163.958659999659, 62.2739300405841, 62.2695850000478, 
    62.2648999998352, 62.2581699997457, 62.2527699998357, 62.2473699999256, 
    62.2422650001297, 62.2371600003338, 62.2304550000483, 62.2237499997627, 
    62.2154199996334, 62.2137999995704, 62.2114799996024, 62.2112999998452, 
    62.211120000088, 62.2094475000771, 62.2077750000662, 62.2061025000554, 
    62.2044300000445), .Dim = c(19L, 2L), class = c("XY", "LINESTRING", 
    "sfg")), structure(c(-162.435819999793, -162.430109999976, 
    -162.427880000112, -162.425650000247, -162.425349999903, 
    -162.42719999978, -162.428489999914, -162.42620000043, -162.421750000262, 
    -162.417300000095, -162.41132999998, -162.405359999864, -162.399389999749, 
    -162.390509999887, -162.381630000024, -162.37460000004, -162.367570000055, 
    -162.363089999854, -162.358609999652, -162.352496666488, 
    -162.346383333324, -162.34027000016, -162.333680000232, -162.329452505243, 
    61.9670499997578, 61.9634399996563, 61.9593549996835, 61.9552699997108, 
    61.9520299995849, 61.9504499995679, 61.9418499995779, 61.9350699998806, 
    61.9318899998237, 61.9287099997667, 61.9272066666544, 61.9257033335421, 
    61.9242000004297, 61.9244350000253, 61.9246699996208, 61.9266249998439, 
    61.928580000067, 61.9307850001277, 61.9329900001884, 61.9375400001709, 
    61.9420900001534, 61.946640000136, 61.9500199999731, 61.9521882750156
    ), .Dim = c(24L, 2L), class = c("XY", "LINESTRING", "sfg"
    )), structure(c(-162.329452505243, -162.327090000303, -162.321173333582, 
    -162.315256666862, -162.309340000141, -162.300583333454, 
    -162.291826666767, -162.28307000008, -162.275282500012, -162.267494999944, 
    -162.259707499876, -162.251919999808, -162.244499999825, 
    -162.237079999842, -162.230319999868, -162.223559999894, 
    -162.21679999992, -162.209842499906, -162.202884999893, -162.195927499879, 
    -162.188969999866, 61.9521882750156, 61.9533999998102, 61.9556199999631, 
    61.957840000116, 61.9600600002689, 61.9619966669127, 61.9639333335566, 
    61.9658700002004, 61.9665250000538, 61.9671799999071, 61.9678349997605, 
    61.9684899996138, 61.9682049999608, 61.9679200003078, 61.9667066667793, 
    61.9654933332508, 61.9642799997222, 61.9621174998228, 61.9599549999234, 
    61.9577925000239, 61.9556300001245), .Dim = c(21L, 2L), class = c("XY", 
    "LINESTRING", "sfg"))), class = c("sfc_LINESTRING", "sfc"
), precision = 0, bbox = structure(c(-164.005879999962, 61.9242000004297, 
-162.188969999866, 62.4387199999281), .Names = c("xmin", "ymin", 
"xmax", "ymax"), class = "bbox"), crs = structure(list(epsg = 4326L, 
    proj4string = "+proj=longlat +datum=WGS84 +no_defs"), .Names = c("epsg", 
"proj4string"), class = "crs"), n_empty = 0L)), .Names = c("id", 
"segnum", "geometry"), row.names = c(NA, 5L), class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_), .Names = c("id", "segnum"), .Label = c("constant", 
"aggregate", "identity"), class = "factor"))

Answer 1

geosphere::dist2Line returns the line ID of a SpatialLines* object.

I have shown that here for the same (or similar) problem as yours:

https://stackoverflow.com/questions/47675571/r-spatial-join-between-spatialpoints-gps-coordinates-and-spatiallinesdatafra

It also shows a nicer approach to providing example data than the output of dput

Answer 2

You can split the sf points object in parts and calculate the distances of every part to lines. I suggest use the library parallel to speed up the process. Also, I think that calculating the distances to lines in small parts use less RAM memory.

My answer is partially based in this post Processing vector to raster faster with R. If you are using Mac OS or Windows OS check this post Efficient spatial joining for large dataset in R to build/configure your cluster.

Try the reproducible and commented example below. At the end you have the nearest line segment ID (segnum) to each point (system.index).

# Load libraries ----------------------------------------------------------

library(sf)
library(geosphere)

# Create sample data ------------------------------------------------------

# 500000 random  points
extentPolygon <- st_as_sfc(st_bbox(lines)) # extent polygon of lines sf object
randomPoints <- st_sample(extentPolygon, size = 500000, type = "random") # randomPoints in extent Polygon
pixels <- st_sf(system.index = 1:500000, geom = randomPoints, crs = 4326) # add data to randomPoints sfc object

# Plot lines + sample points
plot(lines[2], axes = TRUE, graticule = TRUE, xlab = "Longitude", ylab = "Latitude", lwd = 4)
plot(pixels[1:10000,], pch = 19, col = "black", cex = 0.2, add = TRUE) # plot a sample of 10000 points (plotting all the points fill the map)

Without Cluster

# Simple thread dist2Line function ----------------------------------------

# Very time consuming: I tried it with a sample of pixels

points.sp <- as_Spatial(st_geometry(pixels[1:10000,]))
points.sp$system.index <- pixels$system.index[1:10000]
lines.sp <- as_Spatial(st_geometry(lines), IDs = as.character(lines$segnum)) 

system.time(dist <- geosphere::dist2Line(p = points.sp, line = lines.sp))

# user  system elapsed 
# 91.728   0.012  91.752 

# Convert dist to data.frame
dist.df <- as.data.frame(dist)

# Add points ID column to dist.df
dist.df$system.index <- pixels$system.index[1:1000]
colnames(dist.df) <- c("distance", "lon", "lat", "segnum", "system.index")

Using Cluster

# Make cluster ------------------------------------------------------------

# Load 'parallel' package for support Parallel computation in R
library(parallel)

# Calculate the number of cores
no_cores <- detectCores() - 1

# Split features in n parts
n <- 100 # split pixels by 100 points
parts <- split(1:nrow(pixels), cut(1:nrow(pixels), n))

lines.sp <- as_Spatial(st_geometry(lines), IDs = as.character(lines$segnum)) 

# Initiate cluster (after loading all the necessary object to R environment: BRA_adm2, parts, r.raster, n)
cl <- makeCluster(no_cores, type = "FORK")
print(cl)    

# Multithread dist2Line function ------------------------------------------

# Parallelize dist2Line function
system.time(distParts <- parLapply(cl = cl, 
                                   X = 1:n, 
                                   fun = function(x) {
                                     points.sp <- as_Spatial(st_geometry(pixels[parts[[x]],]))
                                     points.sp$system.index <- pixels$system.index[parts[[x]]]
                                     dist <- geosphere::dist2Line(p = points.sp, line = lines.sp)
                                     # Convert dist to data.frame
                                     dist.df <- as.data.frame(dist)
                                     # Add points ID column to dist.df
                                     dist.df$system.index <- pixels$system.index[parts[[x]]]
                                     colnames(dist.df) <- c("distance", "lon", "lat", "segnum", "system.index")
                                     gc(verbose = FALSE) # free memory
                                     return(dist.df)
                                   }))

# user   system  elapsed 
# 1.869    2.067 1314.271

# Finish
stopCluster(cl)

# Bind rows
distBind <- do.call("rbind", distParts)