Statisfaction

Create maps with maptools R package

Posted in Dataset, R by Julyan Arbel on 13 December 2011

Baptiste Coulmont explains on his blog how to use the R package maptools. It is based on shapefile files, for example the ones offered by the French geography agency IGN (at départements and communes level). Some additional material like roads and railways are provided by the OpenStreetMap project, here. For the above map, you need to dowload and dezip the files departements.shp.zip and ile-de-france.shp.zip. The red dots correspond to points of interest longitude / latitude, here churches stored in a vector eglises (use e.g. this to geolocalise places of interest). Then run this code from Baptiste’s tutorial

library(maptools)
france<-readShapeSpatial("departements.shp",
  proj4string=CRS("+proj=longlat"))
routesidf<-readShapeLines(
  "ile-de-france.shp/roads.shp",
  proj4string=CRS("+proj=longlat")
  )
trainsidf<-readShapeLines(
  "ile-de-france.shp/railways.shp",
  proj4string=CRS("+proj=longlat")
  )
plot(france,xlim=c(2.2,2.4),ylim=c(48.75,48.95),lwd=2)
plot(routesidf[routesidf$type=="secondary",],add=TRUE,lwd=2,col="lightgray")
plot(routesidf[routesidf$type=="primary",],add=TRUE,lwd=2,col="lightgray")
plot(trainsidf[trainsidf$type=="rail",],add=TRUE,lwd=1,col="burlywood3")
points(eglises$lon,eglises$lat,pch=20,col="red")

Created by Pretty R at inside-R.org

About these ads

4 Responses

Subscribe to comments with RSS.

  1. dimitris said, on 5 January 2012 at 10:36

    hi and happy new year,
    i tried to reproduce your example because i have some troubles combining openstreet maps and maptools. For some reason i cannot load the shp files from IGN but my french are really poor to understand if i am making a mistake :)

    i was wondering how do you combine these two? if i export an area from openstreet maps, lets say this (http://imageshack.us/photo/my-images/24/screenshotat20120105101.png/) and plot it with maptools with shapefiles for the world like this
    w<-readShapePoly("world/countries.shp",proj4string=CRS("+proj=longlat")
    plot(w,xlim=c(-5.8,29.1),ylim=c(36.3,53))
    i get the export area from openstreet maps and put it xlim and ylim
    i get this result http://imageshack.us/photo/my-images/842/screenshotat201dsa.jpg/

    which is different. From what i see from your example, you use longlat, i use the same projection but i get different plot areas..which projection (i guess it is a projection issue cause i am newbie in plotting maps) should i use to to plot correctly the exported area from openstreet maps with maptools in R?

    dimitris

  2. Julyan Arbel said, on 5 January 2012 at 15:32

    Hi Dimitris, you make me realize there is the same issue for the map above : ylim bounds are correct, but xlim ones are drawn wider than specified. Only a few minutes in my case, but several degrees for your map. Probably an issue of projection, see here for the use of rgdal package (transformation operations): https://stat.ethz.ch/pipermail/r-sig-geo/2008-April/003369.html and here for a nice description of projection systems http://xkcd.com/977/ :)

    • dimitris said, on 7 January 2012 at 11:02

      hi,
      so i am not the only one with this problem :)

      the link from the mailing list gives a reasonable explanation.
      thanks!

  3. Map | m's R Blog said, on 12 June 2013 at 20:39

    […] maptools (homepage, tutorial) […]


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 53 other followers

%d bloggers like this: