StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

POConvert latitude and longitude coordinates to country name in R
text
Body
copied!<p>I have a list of latitude and longitude coordinates, and wish to find out which country they all reside in.</p> <p>I modified an answer from <a href="https://stackoverflow.com/questions/8751497/latitude-longitude-coordinates-to-state-code-in-r">this question about lat-long to US states</a>, and have a working function, but I run into the problem that the <code>worldHires</code> map (from the <code>mapdata</code> package) is hideously out of date and contains a lot of obsolete countries such as Yugoslavia and the USSR.</p> <p>How would I modify this function to use a more modern package, such as <code>rworldmap</code>? I have only managed to frustrate myself so far...</p> <pre><code>library(sp) library(maps) library(rgeos) library(maptools) # The single argument to this function, points, is a data.frame in which: # - column 1 contains the longitude in degrees # - column 2 contains the latitude in degrees coords2country = function(points) { # prepare a SpatialPolygons object with one poly per country countries = map('worldHires', fill=TRUE, col="transparent", plot=FALSE) names = sapply(strsplit(countries$names, ":"), function(x) x[1]) # clean up polygons that are out of bounds filter = countries$x < -180 & !is.na(countries$x) countries$x[filter] = -180 filter = countries$x > 180 & !is.na(countries$x) countries$x[filter] = 180 countriesSP = map2SpatialPolygons(countries, IDs=ids, proj4string=CRS("+proj=longlat +datum=wgs84")) # convert our list of points to a SpatialPoints object pointsSP = SpatialPoints(points, proj4string=CRS("+proj=longlat +datum=wgs84")) # use 'over' to get indices of the Polygons object containing each point indices = over(pointsSP, countriesSP) # Return the state names of the Polygons object containing each point myNames = sapply(countriesSP@polygons, function(x) x@ID) myNames[indices] } ## ## this works... but it has obsolete countries in it ## # set up some points to test points = data.frame(lon=c(0, 5, 10, 15, 20), lat=c(51.5, 50, 48.5, 47, 44.5)) # plot them on a map map("worldHires", xlim=c(-10, 30), ylim=c(30, 60)) points(points$lon, points$lat, col="red") # get a list of country names coords2country(points) # returns [1] "UK" "Belgium" "Germany" "Austria" "Yugoslavia" # number 5 should probably be in Serbia... </code></pre>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload