Search and draw cities on a map using OpenStreetMap and R

(this tutorial can also by used dummies in geography…)

This tutorial explains how to build a map from a matrix containing city names and corresponding country codes, as well as, eventually a frequency attached to each country. The script

• first find thanks to the nice open projet Open Street Map the coordinates of the corresponding cities. This step uses the R package RJSONIO to parse the result of the query performed on open street map;
• from this query, the map is made, using the R package OpenStreetMap (see my previous post for an example of how to use this package and indications on how to install it).

The data

The data I used for this application were coming from Avner (a colleague who must now forget all he knows about this sh*** Go*** Map…) and are highly confidential! I thus provide you a toy example with a data frame with city names, country codes (all cities located in France) and frequencies.

cities = data.frame(nom=c("Toulouse", "Paris", "Égletons", "Marseille", 
                           "Clermont Ferrand"), pays=rep("FR",5),
                     effectif=c(20,5,15,3,3))
print(cities)

that gives the following data frame:

               nom pays effectif
1         Toulouse   FR       20
2            Paris   FR        5
3         Égletons   FR       15
4        Marseille   FR        3
5 Clermont Ferrand   FR        3

Note that some of the names of the cities are a bit tricky to handle (with accents on letters or spaces in the names).

Find the coordinates of the cities

Using the R package RJSONIO, I made a function which, using city name and country code, first sends a query to OpenStreetMap for locating the city and then downloads the city coordinates (in terms of latitude and longitude). City names are simplified replace white spaces with the HTML code “%20”. If the city is not found, the function returns a double NA for both coordinates.

locateCountry = function(nameCity, codeCountry) {
  cleanCityName = gsub(' ', '%20', nameCity)
  url = paste(
    "http://nominatim.openstreetmap.org/search?city="
    , cleanCityName
    , "&countrycodes="
    , codeCountry
    , "&limit=9&format=json"
    , sep="")
  resOSM = fromJSON(url)
  if(length(resOSM) > 0) {
    return(c(resOSM[[1]]$lon, resOSM[[1]]$lat))
  } else return(rep(NA,2)) 
}

The function is then applied to all my cities:

coord = t(apply(cities, 1, function(aRow) locateCountry(aRow[1], aRow[2])))

to obtain a matrix with 5 rows (the cities) and 2 columns (longitude and latitude, in that order).

Draw the map!

This part uses the R package OpenStreetMap.

First, the map is set with proper boundaries (I focused on France, since I only had French countries): the coordinates of the border must be given with two vectors corresponding to the upper left corner and to the down right corner, respectively. I personally use OpenStreetMap to set these values properly:

frenchMap = openmap(c(51.700,-5.669), c(41,11), type="osm")
plot(frenchMap, raster=TRUE)

Then, the coordinates obtained in the previous section are transformed into coordinates suited to OpenStreetMap:

xy = SpatialPoints(apply(coord,2,as.numeric))
proj4string(xy) = CRS("+proj=longlat +ellps=WGS84")
spCities = spTransform(xy, CRS("+proj=merc +a=6378137
  +b=6378137 +lat_ts=0.0+lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m
  +nadgrids=@null +no_defs"))

and all cities are finally represented with dots having a radius proportional to the square root of their frequency. I reduced the alpha chanel of the color to allow for transparency (which is very usefull if some of your cities are overlapping):

sapply(1:nrow(cities), function(ind) {
  plot(spCities[ind,], add=TRUE, col=heat.colors(20,alpha=0.5), lwd=2, pch=19,
       cex=sqrt(as.numeric(cities[ind,3]))/nrow(cities)*2)})

The result is: