R Spatial Join – Nearest Neighbour Join Between Points and Line

linerspatial-join

I want to do spatial join to an OSM line data and GPS point data both in geographic coordinate system, the result which I want is a road layer(has a unique ID) with the attributes of of road and the point file. I don't want to use buffer: sp::over (I get NA values). I tried using gDistance() it only gives me the distance but does not join the attribute tables. I want to replicate NNJoin of QGIS in R since I can't use batch processing for NNJoin.

I don't think I am satisfied by the answer here https://stackoverflow.com/questions/47675571/r-spatial-join-between-spatialpoints-gps-coordinates-and-spatiallinesdatafra. 

longitude <- c(10.86361, 10.96062, 10.93032, 10.93103, 10.93212)        
latitude <- c(44.53355, 44.63234, 44.63470, 44.63634, 44.64559)
longlat <- data.frame(longitude, latitude)  
my_spdf<- SpatialPoints(longlat, proj4string=crs("+init=epsg:4326"))  
shapefile_roads<- spLines(cbind(longitude+.1, latitude), cbind(longitude-.1, rev(latitude)), cbind(longitude-.1, latitude+1), crs=crs("+init=epsg:4326"))
shapefile_roads <- spTransform(shapefile_roads, CRS("+proj=utm +zone=43 ellps=WGS84"))
my_spdf <- spTransform(my_spdf, CRS("+proj=utm +zone=43 ellps=WGS84"))
result <- over(shapefile_roads, my_spdf)

Best Answer

Here's a solution using library(sf):

library(sf)
library(rgdal)
library(sp)
library(raster)

longitude <- c(10.86361, 10.96062, 10.93032, 10.93103, 10.93212)        
latitude <- c(44.53355, 44.63234, 44.63470, 44.63634, 44.64559)
longlat <- data.frame(longitude, latitude)  
my_spdf<- SpatialPoints(longlat, proj4string=CRS("+init=epsg:4326"))  
shapefile_roads<- spLines(cbind(longitude+.1, latitude), cbind(longitude-.1, rev(latitude)), cbind(longitude-.1, latitude+1), crs=CRS("+init=epsg:4326"))
shapefile_roads <- spTransform(shapefile_roads, CRS("+proj=utm +zone=43 ellps=WGS84"))
my_spdf <- spTransform(my_spdf, CRS("+proj=utm +zone=43 ellps=WGS84"))

my_sfdf = st_sf(id_pt = 1:length(my_spdf), geom = st_as_sfc(my_spdf))
shapefile_roads_sf = st_sf(id_ln = 1:length(shapefile_roads), 
                           geom = st_as_sfc(shapefile_roads))

# get nearest point for each line
shapefile_roads_sf_w_pts = st_join(shapefile_roads_sf, 
                                   my_sfdf, 
                                   join = st_nearest_feature)

shapefile_roads_sf_w_pts
Simple feature collection with 3 features and 2 fields
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: -4358172 ymin: 7326658 xmax: -4226019 ymax: 7428067
epsg (SRID):    32643
proj4string:    +proj=utm +zone=43 +ellps=WGS84 +units=m +no_defs
  id_ln id_pt                           geom
1     1     2 LINESTRING (-4350438 732675...
2     2     1 LINESTRING (-4348820 735297...
3     3     5 LINESTRING (-4243179 742480...

Please, when posting a reproducible example, make sure that it does indeed run and also include the packages that need to be loaded in order to execute the code. We cannot know where functions come from if the execution fails! Also, if you say you want IDs, please make sure that your example includes them!

Related Question