4

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

A list of coordinates a.k.a pointsCoordinates

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

enter image description here

5

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

Colors stand for point ID

  • 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... – Tacticus Jul 12 '16 at 8:41
  • @Tacticus, see edits. – dof1985 Jul 12 '16 at 10:47
  • 1
    If I were able to upvote this futher, I'm absolutely certain I would do it. – Tacticus Jul 12 '16 at 11:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.