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

From: Roger Bivand <Roger.Bivand_at_nhh.no>
Date: Wed 10 Jan 2007 - 19:38:05 GMT


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@nhh.no

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Thu Jan 11 06:44:28 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Wed 10 Jan 2007 - 20:30:29 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.