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

[Tord Snäll]

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.

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 ...



-- 

Tord Snäll
Department of Conservation Biology
Swedish University of Agricultural Sciences (SLU)
P.O. 7002, SE-750 07 Uppsala, Sweden
Office/Mobile/Fax
+46-18-672612/+46-730-891356/+46-18-673537
E-mail: tord.snall at nvb.slu.se
www.nvb.slu.se/staff_tordsnall