[R] map data.frame() data after having linked them to a read.shape() object

[Roger Bivand]

On Wed, 10 Jan 2007, Tord Snäll wrote:

> Dear all,
> I try to first link data in a data.frame() to a (polygon) read.shape() 
> object and then to plot a polygon map showing the data in the 
> data.frame. The first linking is called "link" in ArcView and "relate" 
> in ArcMap. I use the code shown below, though without success.
> 
> Help with this would be greatly appreciated.

Please consider continuing this thread on the R-sig-geo list. Searching 
the archives might also have given you some leads. For now, and apart from 
not using sp classes (see note in R News in 2005), you have 490 polygons 
in the shapefile - probably one duplicate and 489 unique entities in 
STACOUNT4, but only 106 unique stacount of 431 observations in the data 
frame. The plot method for Map objects is deprecated. The classes in the 
sp package use S4, not S3, specifically to help users avoid putting things 
inside objects that break methods.

In maptools, see ?readShapePoly, and ?spCbind to see how to read a
shapefile into an sp object (check the 490/489 issue), and how the
Polygons IDs can be used as a key with the data frame row names to make
this easier to do. Please also consider using FIPS numbers as IDs for
counties; the five digit ssccc ID is fairly standard, and avoids the
problem of repetitive county names across US states.

> 
> Thanks!
> 
> Tord
> 
> require(maptools)
> # Read shape file (one row per county)
> a=read.shape("myshp.shp", dbf.data=TRUE, verbose=TRUE)
> str(a)
>   ..- attr(*, "maxbb")= num [1:4] -100   49    0    0
>   ..- attr(*, "class")= chr "ShapeList"
> $ att.data:'data.frame':       490 obs. of  60 variables:
>   ..$ STATE_FIPS: Factor w/ 12 levels "04","06","08",..: 11 11 11 11 4 5
> 5 5 5 5 ...
> [snip]
>   ..$ STACOUNT4 : Factor w/ 489 levels "ArizonaApache",..: 437 460 451
> 453 147 207 195 198 231 206 ...
>   ..- attr(*, "data_types")= chr [1:60] "C" "C" "C" "C" ...
> - attr(*, "class")= chr "Map"
> 
> 
> # Read case data (one row per case)
> cases = read.table("cases.txt", h=T,)
> str(cases)
> 'data.frame':   431 obs. of  8 variables:
> $ Year    : int  1950 1950 1950 1951 1956 1957 1959 1959 1959 1959 ...
> $ Case    : int  3 1 2 1 1 1 2 4 1 3 ...
> $ stacount: Factor w/ 106 levels "ArizonaApache",..: 1 66 76 66 26 29
> 15 25 30 60 ...
> 
> # table the cases data PER Year, PER County (County = "stacount")
> temp = t(table(cases[,c("Year","stacount")]))
> stacount = dimnames(temp)$stacount
> temp = cbind(stacount, as.data.frame(temp[,1:ncol(temp)],row.names=F))
> str(temp)
> 'data.frame':   106 obs. of  50 variables:
> $ stacount: Factor w/ 106 levels "ArizonaApache",..: 1 2 3 4 5 6 7 8 9
> 10 ...
> $ 1950    : int  1 0 0 0 0 0 0 0 0 0 ...
> [snip]
> $ 2005    : int  0 0 0 0 0 0 0 0 0 0 ...
> 
> # Pick out a temporary attribute data.frame
> datfr = a$att.data
> 
> # Merge the temporaty data frame with tabled "cases"
> for(i in 2:ncol(temp)){
>      datfr = merge(datfr, temp[,c(1,i)], by.x="STACOUNT4",
> by.y="stacount", all.x=T, all.y=F)
> }
> 
> #Replace NAs with 0:
> for(i in 61:109){
>      datfr[,i] = ifelse(is.na(datfr[,i])==T,0,datfr[,i])
> }
> 
> str(a$att.data)
> 'data.frame':   490 obs. of  60 variables:
> $ NAME      : Factor w/ 416 levels "Ada","Adams",..: 120 352 265 277 33
> 210 122 135 372 209 ...
> [snip]
> $ STACOUNT4 : Factor w/ 489 levels "ArizonaApache",..: 437 460 451 453
> 147 207 195 198 231 206 ...
> - attr(*, "data_types")= chr  "C" "C" "C" "C" ...
> # Note that the above data is of "attribute type"
> 
> str(datfr)
> 'data.frame':   490 obs. of  109 variables:
> $ STACOUNT4 : Factor w/ 489 levels "ArizonaApache",..: 1 2 3 4 5 6 7 8
> 9 10 ...
> [snip]
> $ 1951      : num  0 0 0 0 0 0 0 0 0 0 ...
> [snip]
> $ 2005      : num  0 0 0 0 0 0 0 0 0 0 ...
> # Note that at the end of this, data type is not described - it is a 
> "simple" data frame
> 
> # bind data together:
> #Alternative 1:
> a$att.data = cbind(a$att.data, datfr[,61:109])
> # Other alternatives:
> test = matrix(ncol=49)
> a$att.data[,61:109] = test
> a$att.data[,61:109] = datfr[,61:109]
> 
> # plot:
> plot(a, auxvar=a$att.data[,61], xlim=c(-125,-99),ylim=c(28,52), xlab="",
> ylab="", frame.plot=F,axes=F)
> There were 50 or more warnings (use warnings() to see the first 50)
> warnings()
> 49: "axes" is not a graphical parameter in: polygon(xy$x, xy$y, 
> col,border, lty, ...)
> 50: "frame.plot" is not a graphical parameter in: polygon(xy$x, 
> xy$y,col, border, lty, ...)
> 
> # The a$att.data type has changed to becoming a typical data.frame - 
> "attr" is not mentioned:
> str(a$att.data)
> [snip]
> $ 2003      : num  0 0 0 0 0 0 0 0 0 0 ...
> $ 2004      : num  0 0 0 0 0 0 0 0 0 0 ...
> $ 2005      : num  0 0 0 0 0 0 0 0 0 0 ...
> 
> 
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no