[GIS] Road snapping in R – snapPointsToLine / Projection Problem

Question

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)
pointsSpatialCorrected = snapPointsToLines(pointsSpatial, roadsBavaria)

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?

---

3. Additional information

The resulting value of spatialPointsCorrected as displayed by View()

Answer 1

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


bavariaUTM <- spTransform(readOGR(".", "roads"), crs) # get road data and project to UTM.

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).