10

I am trying to use the function fortify in R to plot the map. But I am always getting an error while processing the data. I read the shape file using the following code:

shape <- readShapeSpatial("Lond_City/lond_city.shp")

The shapefile can be found on the following dropbox location: https://www.dropbox.com/sh/d4w6saailydtu1r/5oIa56xV6S

Then I tried to use the following command :

shape1 <- fortify(shape,region="data")

In the above code, what should I put in Place of "data". This is the place which seems confusing. I can open the data or shapefile in R. However, when I run the fortify line I get the following error:

Error in '[.data.frame'(attr, , region) : undefined columns selected

Can anybody tell me how to properly fortify a shape file with example ? The reason I want to use fortify is that I want to combine the shape file with the data points.

Thank you so much.

Jdbaba

Jd Baba
  • 5,948
  • 18
  • 62
  • 96
  • It looks like you're using something from the Ordnance Survey? If you can make a link to the shapefile so that others can download it, or explain **exactly** where you got the data from, then I think we can probably help. – SlowLearner Jan 17 '13 at 02:35
  • Slowlearner, thank you for your reply. I have added the link of the shapefile in the question. Please have a look. thanks. – Jd Baba Jan 17 '13 at 02:45
  • Thank you for that. You will nearly always get a faster and better response on StackOverflow if you provide data as well as code. – SlowLearner Jan 17 '13 at 02:59

2 Answers2

12

Your problem is that you need to give fortify a column that actually exists. If you type str(lon.shape) from my example code below, you will see that there is no region column:

> str(lon.shape)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 1 obs. of  7 variables:
  .. ..$ ons_label: chr "00AA"
  .. ..$ name     : chr "City of London"
  .. ..$ label    : chr "02AA"
  .. ..$ X_max    : num 533843
  .. ..$ y_max    : num 182198
  .. ..$ X_min    : num 0
  .. ..$ y_min    : num 0
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 532464 181220
  .. .. .. .. .. .. ..@ area   : num 3151465
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:771, 1:2] 531027 531029 531036 531074 531107 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 532464 181220
  .. .. .. ..@ ID       : chr "0"
  .. .. .. ..@ area     : num 3151465
  ..@ plotOrder  : int 1
  ..@ bbox       : num [1:2, 1:2] 530967 180404 533843 182198
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. ..@ projargs: chr "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs" 

Instead, try using name, like this:

lon.df <- fortify(lon.shape, region = "name")

Which produces this data frame:

> head(lon.df)
      long      lat order  hole piece            group             id
1 531026.9 181611.1     1 FALSE     1 City of London.1 City of London
2 531028.5 181611.2     2 FALSE     1 City of London.1 City of London
3 531036.1 181611.5     3 FALSE     1 City of London.1 City of London
4 531074.0 181610.3     4 FALSE     1 City of London.1 City of London
5 531107.0 181609.3     5 FALSE     1 City of London.1 City of London
6 531117.1 181608.9     6 FALSE     1 City of London.1 City of London

Here's one way to go about it, start to finish:

library(rgdal)
library(ggplot2)
library(rgeos)

shape.dir <- "c:/test/london" # use your directory name here

lon.shape <- readOGR(shape.dir, layer = "lond_city")
lon.df <- fortify(lon.shape, region = "name")

ggplot(lon.df, aes(x = long, y = lat, group = group)) +
    geom_polygon(colour = "black", fill = "grey80", size = 1) +
    theme()

...resulting in the following:

screenshot

SlowLearner
  • 7,907
  • 11
  • 49
  • 80
  • SlowLearner, Thank you so much for your wonderful solution. But it seems like I cannot run because I get the following error "isTRUE(gpclibPermitStatus()) is not TRUE. Do you know how to fix this ? I tried to install "gpclibPermit" but it says this package is not available for version 2.15.2. Any suggestions ? – Jd Baba Jan 17 '13 at 03:23
  • 1
    @Jdbaba This is because SlowLearner is using ```readOGR``` rather than ```readShapeSpatial```. You can fix this by running gpclibPermit() but be careful as this is forbidden if you are using the software for a commercial purpose. You may prefer to avoid this issue by adopting the approach in my answer. Alternatively you can install the ```rgeos``` package which has no commercial restriction and should also fix the problem. – orizon Jan 17 '13 at 03:38
  • 1
    @Jdbaba `readOGR` preserves the proj4 information so I generally prefer it. I have edited my answer a little to include `gpclib` and the `gpclibPermit()` call. – SlowLearner Jan 17 '13 at 03:48
  • @SlowLearner If you include ```library(rgeos)``` you do not need to run ```gpclibPermit()``` as well, since ```rgeos``` is an alternative to gpclib. – orizon Jan 17 '13 at 03:54
  • @orizon Ah, thanks for that, guess that's why in another piece of code I have there **is** an unexplained call to `rgeos`! I get confused as to which libraries supply what: `sp`, `map`, `maptools`, `gpclib`, `rgdal`, `rgeos` all get mixed up together in my head. Will edit. – SlowLearner Jan 17 '13 at 03:57
  • Thanks guys for clearing things up. I have another task to do. I need to add the points and contour to the plot and clip the contour by the shapefile. Is that possible in R ? Thanks. – Jd Baba Jan 17 '13 at 04:04
  • @Jdbaba On StackOverflow we like to do things "one question at a time" so please make this a separate question and include the data (or fake data that resembles the real data). – SlowLearner Jan 17 '13 at 04:06
  • @Jdbaba Yes it is possible, search for answers on http://gis.stackexchange.com/ as well as StackOverflow. – SlowLearner Jan 17 '13 at 04:12
2

You should put label in place of data.

require(maptools); require(ggplot2)
shape <- readShapeSpatial("Lond_City/lond_city.shp")
shape1 <- fortify(shape,region="label")

To find this out I inspected the data element of shape:

> shape@data
  ons_label           name label    X_max    y_max X_min y_min
0      00AA City of London  02AA 533842.7 182198.4     0     0

This shows ons_label or name would also have worked and may be more suitable for your purposes.

Note that your shape file is slightly unusual because there seems to be only one group.

orizon
  • 3,159
  • 3
  • 25
  • 30