combine_maptools.R
In maptools: Tools for Handling Spatial Objects

### R code from vignette source 'combine_maptools.Rnw'

###################################################
### code chunk number 1: combine_maptools.Rnw:44-47
###################################################
owidth <- getOption("width")
options("width"=90)
.PngNo <- 0


###################################################
### code chunk number 2: afig (eval = FALSE)
###################################################
## .PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="")
## pdf(file=file, width = 6.5, height = 3.5, pointsize = 12, bg = "white")
## opar <- par(mar=c(3,3,1,1)+0.1)


###################################################
### code chunk number 3: afigl (eval = FALSE)
###################################################
## .PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="")
## pdf(file=file, width = 6.5, height = 3.5, pointsize = 12, bg = "white")


###################################################
### code chunk number 4: bfigl (eval = FALSE)
###################################################
## .PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="")
## pdf(file=file, width = 6.5, height = 5, pointsize = 12, bg = "white")


###################################################
### code chunk number 5: bfig (eval = FALSE)
###################################################
## .PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="")
## pdf(file=file, width = 6.5, height = 5, pointsize = 12, bg = "white")
## opar <- par(mar=c(3,3,1,1)+0.1)


###################################################
### code chunk number 6: zfig (eval = FALSE)
###################################################
## par(opar)
## dev.null <- dev.off()
## cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")


###################################################
### code chunk number 7: zfigl (eval = FALSE)
###################################################
## dev.null <- dev.off()
## cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")


###################################################
### code chunk number 8: combine_maptools.Rnw:101-111
###################################################
owd <- getwd()
setwd(system.file("shapes", package="maptools"))
library(maptools)
nc90 <- readShapeSpatial("co37_d90")
proj4string(nc90) <- CRS("+proj=longlat +datum=NAD27")
sc90 <- readShapeSpatial("co45_d90")
proj4string(sc90) <- CRS("+proj=longlat +datum=NAD27")
va90 <- readShapeSpatial("co51_d90")
proj4string(va90) <- CRS("+proj=longlat +datum=NAD27")
setwd(owd)


###################################################
### code chunk number 9: combine_maptools.Rnw:117-124
###################################################
.PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="")
pdf(file=file, width = 6.5, height = 5, pointsize = 12, bg = "white")
opar <- par(mar=c(3,3,1,1)+0.1)
oopar <- par(mar=c(3,2,1,1)+0.1)
plot(va90, xlim=c(-85,-75), ylim=c(32,40), axes=TRUE, border="grey10")
plot(nc90, add=TRUE, border="grey40")
plot(sc90, add=TRUE, border="grey70")
par(oopar)
par(opar)
dev.null <- dev.off()
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")


###################################################
### code chunk number 10: combine_maptools.Rnw:148-149
###################################################
library(maptools)


###################################################
### code chunk number 11: combine_maptools.Rnw:151-155
###################################################
names(sc90)
sc90a <- spChFIDs(sc90, paste(sc90$ST, sc90$CO, sep=""))
sc90a <- sc90a[,-(1:4)]
names(sc90a)


###################################################
### code chunk number 12: combine_maptools.Rnw:157-158
###################################################
proj4string(sc90a) <- CRS(proj4string(sc90a))


###################################################
### code chunk number 13: combine_maptools.Rnw:169-170
###################################################
names(nc90)


###################################################
### code chunk number 14: combine_maptools.Rnw:172-173 (eval = FALSE)
###################################################
## nc90a <- spChFIDs(nc90, paste(nc90$ST, nc90$CO, sep=""))


###################################################
### code chunk number 15: combine_maptools.Rnw:175-176
###################################################
try1 <- try(spChFIDs(nc90, paste(nc90$ST, nc90$CO, sep="")), silent=TRUE)


###################################################
### code chunk number 16: combine_maptools.Rnw:178-179
###################################################
cat(try1)


###################################################
### code chunk number 17: combine_maptools.Rnw:189-190
###################################################
table(table(paste(nc90$ST, nc90$CO, sep="")))


###################################################
### code chunk number 18: combine_maptools.Rnw:216-219
###################################################
if (rgeosStatus()) {
nc90a <- unionSpatialPolygons(nc90, IDs=paste(nc90$ST, nc90$CO, sep=""))
}


###################################################
### code chunk number 19: combine_maptools.Rnw:229-234
###################################################
if (rgeosStatus()) {
nc90_df <- as(nc90, "data.frame")[!duplicated(nc90$CO),-(1:4)]
row.names(nc90_df) <- paste(nc90_df$ST, nc90_df$CO, sep="")
nc90b <- SpatialPolygonsDataFrame(nc90a, nc90_df)
}


###################################################
### code chunk number 20: combine_maptools.Rnw:254-263
###################################################
if (rgeosStatus()) {
va90a <- spChFIDs(va90, paste(va90$ST, va90$CO, sep=""))
va90a <- va90a[,-(1:4)]
va90_pl <- slot(va90a, "polygons")
va90_pla <- lapply(va90_pl, checkPolygonsHoles)
p4sva <- CRS(proj4string(va90a))
vaSP <- SpatialPolygons(va90_pla, proj4string=p4sva)
va90b <- SpatialPolygonsDataFrame(vaSP, data=as(va90a, "data.frame"))
}


###################################################
### code chunk number 21: combine_maptools.Rnw:303-309
###################################################
if (rgeosStatus()) {
nc_sc_va90 <- spRbind(spRbind(nc90b, sc90a), va90b)
FIPS <- row.names(nc_sc_va90)
str(FIPS)
length(slot(nc_sc_va90, "polygons"))
}


###################################################
### code chunk number 22: combine_maptools.Rnw:338-340
###################################################
t1 <- read.fwf(system.file("share/90mfips.txt", package="maptools"), skip=21,
 widths=c(4,4,4,4,2,6,2,3,3,1,7,5,3,51), colClasses = "character")


###################################################
### code chunk number 23: combine_maptools.Rnw:342-348
###################################################
t2 <- t1[1:2004,c(1,7,8,14)]
t3 <- t2[complete.cases(t2),]
cnty1 <- t3[t3$V7 != "  ",]
ma1 <- t3[t3$V7 == "  ",c(1,4)]
cnty2 <- cnty1[which(!is.na(match(cnty1$V7, c("37", "45", "51")))),]
cnty2$FIPS <- paste(cnty2$V7, cnty2$V8, sep="")


###################################################
### code chunk number 24: combine_maptools.Rnw:368-375
###################################################
if (rgeosStatus()) {
MA_FIPS <- cnty2$V1[match(FIPS, cnty2$FIPS)]
MA <- ma1$V14[match(MA_FIPS, ma1$V1)]
MA_df <- data.frame(MA_FIPS=MA_FIPS, MA=MA, row.names=FIPS)
nc_sc_va90a <- spCbind(nc_sc_va90, MA_df)
ncscva_MA <- unionSpatialPolygons(nc_sc_va90a, nc_sc_va90a$MA_FIPS)
}


###################################################
### code chunk number 25: combine_maptools.Rnw:381-392
###################################################
.PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="")
pdf(file=file, width = 6.5, height = 5, pointsize = 12, bg = "white")
opar <- par(mar=c(3,3,1,1)+0.1)
if (rgeosStatus()) {
oopar <- par(mar=c(3,2,1,1)+0.1)
plot(nc_sc_va90, border="grey", axes=TRUE)
plot(ncscva_MA, lwd=2, add=TRUE)
text(coordinates(ncscva_MA), labels=row.names(ncscva_MA), cex=0.6)
par(oopar)
} else {
plot(1)
}
par(opar)
dev.null <- dev.off()
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")


###################################################
### code chunk number 26: combine_maptools.Rnw:410-417
###################################################
if (rgeosStatus()) {
np <- sapply(slot(ncscva_MA, "polygons"), function(x) length(slot(x, "Polygons")))
table(np)
MA_fips <- row.names(ncscva_MA)
MA_name <- ma1$V14[match(MA_fips, ma1$V1)]
data.frame(MA_fips, MA_name)[np > 1,]
}


###################################################
### code chunk number 27: combine_maptools.Rnw:425-426
###################################################
options("width"=owidth)