Incorrect NA return when converting Lat/Long Coordinates to location in R

Question

I am trying to use a modified version of the R code found in the following link:

Latitude Longitude Coordinates to State Code in R

To test the code, I created the following formal arguments:

mapping = "state"
pointsDF = data.frame(x = c(-88.04607, -83.03579), y = c(42.06907, 42.32983))
latlong2state(pointsDF, mapping)

The code returned the following:

[1] "Illinois" NA

The first coordinate set returns a correct answer, i.e. "Illinois". However, when I input the 2nd coordinate set (i.e. -83.03579, 42.32983) into an online converter, I get the following:

Downtown, Detroit, MI, USA

(http://www.latlong.net/Show-Latitude-Longitude.html)

Running the code again but changing the second coordinate from 42.32983 to 43.33 puts the point in the state of Michigan.

When using the "world" map as my formal argument for the "mapping" variable, the code returns "USA". I have been struggling for days to figure this out and have had no luck. I have played around with SpatialPointDataFrames, various projections, and looked into the state polygon objects themselves. I am using R version 3.3.1 on a Windows 7 system. I think the data point in question may be falling on a border line. In which case, I think an "NA" would be expected. The code I used is below.

Code Used:

library(sp)

library(maps)  
library(maptools)  
library(rgdal)

latlong2state = function(pointsDF, mapping) {

        local.map = map(database = mapping, fill = TRUE, col = "transparent", plot = FALSE) 
        IDs = sapply(strsplit(local.map$names, ":"), function(x) x[1])
        maps_sp = map2SpatialPolygons(map = local.map, ID = IDs, 
                                      proj4string = CRS("+proj=longlat +datum=WGS84"))                        
        pointsSP = SpatialPoints(pointsDF, 
                                 proj4string = CRS("+proj=longlat +datum=WGS84"))
        indices = over(x = pointsSP, y = maps_sp)
        mapNames = sapply(maps_sp@polygons, function(x) {x@ID})
        mapNames[indices]
}

I am only two months in to learning R and love the language thus far. This has been the first time I could not find an answer. I would really appreciate an help provided on the matter!!!

Answer 1

Firstly, the issue is not due to the point lying on a border. In fact, over() would not return NA for a point on a border, but rather "if a point falls in multiple polygons, the last polygon is recorded."

NA denotes a point that does not fall in a polygon. We can zoom in on your map to see this is the case

plot(local.map,  xlim = c(-83.2, -82.8), ylim=c(42.2,42.6), type="l")
polygon(local.map, col="grey60")
points(local.map)
points(pointsDF[2,], col="red")

The point falls outside the contiguous USA in Canada, according to the polygons provided by maps::map(). Why would this be the case when other maps, as you say, locate this point on the USA side of the border? I do not think this is a projection issue, because we are using the same WGS84 geographic coordinates for the polygons and the points. It seems, therefore, that the polygons themselves that are provided by maps::map() may be wrong.

We can check this by comparing to polygons from another source. I downloaded the US census departments highest resolution state boundaries from http://www2.census.gov/geo/tiger/GENZ2015/shp/cb_2015_us_state_500k.zip. Then,

shp.path <- "C:/Users/xxx/Downloads/cb_2015_us_state_500k/cb_2015_us_state_500k.shp"
states <- readOGR(path.expand(shp.path), "cb_2015_us_state_500k")
plot(states,  xlim = c(-83.2, -82.8), ylim=c(42.2,42.6))
points(pointsDF[2,], col="red")

gets us this map in which we see that the point is inside the US boundary:

The solution I recommend therefore, is to use these better resolution, more reliable boundary polygons, particularly if you are interested to accurately resolve points close to borders.