# Road snapping in R - snapPointsToLine / Projection Problem

1. Introduction (you may skip this if you like)

In order to implement a road snapping feature in R for recorded positioning data, I first searched the internet for possible solutions to the problem:

While it seems this can be solved with fuzzyMM (using fuzzy logic), the project using that approach is unfortunately not active any more and, thus, the documentation does not include the kind of information that would be required to get started with the package itself.

Another solution could be using online geospatial APIs, e.g. Google's. However, the current database already includes way more datasets than what's included in the free package.

At last, I found `maptool`'s function `snapPointsToLines` (short introduction about how the function works can be found here) and attempted to snap a couple of coordinates to lines provided by a shapefile.

Note: Because the GPS-data was recorded on cars mostly driving on major roads in Germany, I downloaded the respective shapefile (of Bavaria, i.e. Bayern) here (`.shp.zip`).

2. The specific problem

Road snapping with `snapPointsToLines()` generates quite confusing values with this code. The following picture shows a list of recorded longitude / latitude values (see `pointsCoord`).

``````roadsBavaria = readOGR("~/roads.shp", "roads")
pointsCoord  = getValuesFromPicture()

pointsSpatial          = SpatialPoints(pointsCoord)
``````

Although this reduced test case sets the value `pointsSpatialCorrected` (see result at the end), it fails (code after that is not executed) with

``````Warning messages:
1: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") :
Spatial object 1 is not projected; GEOS expects planar coordinates
2: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") :
Spatial object 2 is not projected; GEOS expects planar coordinates
``````

Note that this error occurs even when I apply the shapefile's projection on the `spatialPoints` manually via

``````pointsSpatial = SpatialPoints(pointsCoord, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
``````

What am I doing wrong here?

The resulting value of `spatialPointsCorrected` as displayed by `View()`

The warning message is quite clear here:

``````1: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") :
Spatial object 1 is not projected; GEOS expects planar coordinates
``````

It requires planar (Cartesian) coordinates, i.e. in meters, miles, etc. You however use polar reference system (WGS-1984).

You should re-project your SpatialPointDataFrame to a Cartesian projection of-your-choice (UTM is one option) and run the function again - it should fix this error. See this thread for additional info on re-projection to UTM: converting-latitude-and-longitude-points-to-utm

EDIT: Reprojecting to UTM zone = 32N, code and results

See the following code and image. Re-projecting correctly fixes the problem. The code using the data provided by the OP.

``````library(maptools)
library(rgdal)
library(sp)

crs <- CRS("+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")     # UTM zone = 32 N
wgs84 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")  # To create Points

coords <- cbind("lon" = c(9.967067, 9.967218,
9.967507, 9.967502,
9.967490, 9.967453),
"lat" = c(50.20421, 50.20419,
50.20411, 50.20410,
50.20410, 50.20410))

pts <- spTransform(SpatialPointsDataFrame(coords = coords,
proj4string = wgs84,
data = data.frame("id" = 1:6)), crs) # create points and transform

pointsSpatialCorrected <- snapPointsToLines(pts, bavariaUTM)      # snap to the nearest road
``````

Here are the original points coordinates in meters:

``````pts@coords
lon     lat
[1,] 569012.3 5561784
[2,] 569023.1 5561782
[3,] 569043.8 5561773
[4,] 569043.5 5561772
[5,] 569042.6 5561772
[6,] 569040.0 5561772
``````

and the coordinates of the snapped points:

``````pointsSpatialCorrected@coords
X       Y
[1,] 569072.5 5561820
[2,] 569076.4 5561814
[3,] 569090.1 5561759
[4,] 569089.8 5561757
[5,] 569089.7 5561757
[6,] 569089.5 5561756
``````

I also exported shp files of these two objects using: `writeOGR(pointsSpatialCorrected, dsn = ".", layer = "correctedPTS.shp", driver = "ESRI Shapefile")` and produced this image in QGIS (It can be plotted in R as well using `plot`, but the roads takes too much time).

• Umm yesterday I didn't really take a look at the values returned, but as it turns out the values (now in UTM; converted back to WGS-84 and check via Bing / Google Maps) are even more inaccurate than what I had without reprojecting both shape file and points to UTM. I will try to update my question accordingly... Jul 12, 2016 at 8:41
• @Tacticus, see edits. Jul 12, 2016 at 10:47
• If I were able to upvote this futher, I'm absolutely certain I would do it. Jul 12, 2016 at 11:40