Tord Snäll
2007-Jan-10 16:49 UTC
[R] map data.frame() data after having linked them to a read.shape() object
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
Roger Bivand
2007-Jan-10 19:38 UTC
[R] map data.frame() data after having linked them to a read.shape() object
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