#######################################################################
# Converts a shapefile to KML polygons for display in Google Maps or
# Google Earth. Adds information on region names and posterior means
# of logistic parameters to a clickable label for each region.
#
# Carl Schmertmann
# created 13 May 10
# altered 14 May 10
#######################################################################
memory.limit(4000)
graphics.off()
library(spdep)
myDir = "."
myShapeFile = "brazil_microregions_1960_2000.shp"
setwd(myDir)
### read the shapefile data for mapping
map = readShapePoly(myShapeFile)
### set up Google Earth-style map for KML information
map@proj4string = CRS("+proj=longlat +ellps=WGS84")
GoogleMap = GE_SpatialGrid(map)
#######################################################################
# INPUT DATA: regional data and summaries of posterior distributions
#######################################################################
## regional data
D = read.csv("brazil_logistic_data.csv")
## load matrix Q containing a summary of posteriors
## Q is an output variable from the WinBUGS program
summary.file = "brazil_logistic_posterior_summary.11-May-10-2200.gz"
load(summary.file)
#############################################
# construct summary labels for clickable map
#############################################
nice = function(x,nplaces=1) { format(round(x,nplaces),nsmall=nplaces) }
label = function(x) { paste("n1980 = ",x[1],"
",x[2]," -",x[3]," ",x[4]," +",x[5],sep="") }
tmp = matrix(NA, nrow(D), 5,
dimnames=list(NULL, c("n80","Fpre","Fpost","start","duration")))
polyname = apply(D[,c("index","name")], 1, paste, collapse=" " )
tmp[,"n80"] = D[,"n80"]
tmp[,"Fpre"] = nice( Q["mean",grep("Fpre", dimnames(Q)[[2]])] , 1 )
tmp[,"Fpost"] = nice( Q["mean",grep("Fpost", dimnames(Q)[[2]])] , 1 )
tmp[,"start"] = floor( 1960+ 10* Q["mean",grep("start", dimnames(Q)[[2]])] )
tmp[,"duration"] = floor( 10* Q["mean",grep("duration", dimnames(Q)[[2]])] )
# lab is a vector of labels with post. means - eg "5.6-2.0 1962 +31"
lab = apply(tmp,1,label)
rm(tmp)
######################################
## colors from PDR maps
hue = c(rgb(138, 66, 38,191,maxColorValue=255),
rgb( 0,163,178,191,maxColorValue=255),
rgb(171,193,121,191,maxColorValue=255),
rgb(248,220,143,191,maxColorValue=255)
)
#######################################
# list of variables, categories, etc
# for mapping
#
# Technical note:
# avoid > or < signs in the legends, because they
# will mess up the KML
#######################################
vars = list(
list(name ="start",
brks = (c(-Inf,1960,1970,1980,Inf)-1960) / 10,
legends = c("before 1960","1960-69","1970-79","after 1980"),
style = 1:4,
pretty.name = "Transition Start Date"),
list(name ="Fpre",
brks = c(-Inf, 5,6,7, Inf),
legends = c("below 5", "5.00-5.99", "6.0-6.99", "above 7"),
style = 1:4,
pretty.name = "Initial Fertility Level Fpre"),
list(name ="Fpost",
brks = c(-Inf, 2, 2.25, 2.5, Inf),
legends = c("below 2", "2.00-2.24","2.25-2.50","above 2.50"),
style = 1:4,
pretty.name = "Final Fertility Level Fpost"),
list(name ="duration",
brks = c(-Inf,20,25,30, Inf)/10 ,
legends = c("less than 20 yrs","20-24 yrs","25-29 yrs","more than 30 years"),
style = 4:1,
pretty.name = "Transition Duration")
)
############################################
# modified kmlPolygon function from maptools
# library
#
# adds altitudes, eliminates separate style
# for each polygon, etc.
############################################
my.kmlPolygon = function (obj = NULL, kmlfile = NULL, name = "R Polygon",
description = "", col = NULL, visibility = 1, lwd = 1, border = 1,
kmlname = "", kmldescription = "", styleNum=1)
{
if (is.null(obj))
return(list(header = c("",
"",
"", paste("", kmlname, "",
sep = ""), paste("", sep = "")), footer = c("",
"")))
if (class(obj) != "Polygons" && class(obj) != "SpatialPolygonsDataFrame")
stop("obj must be of class 'Polygons' or 'SpatialPolygonsDataFrame' [package 'sp']")
if (class(obj) == "SpatialPolygonsDataFrame") {
if (length(obj@polygons) > 1)
warning(paste("Only the first Polygons object with the ID '",
obj@polygons[[1]]@ID, "' is taken from 'obj'",
sep = ""))
obj <- obj@polygons[[1]]
}
kml <- ""
kmlHeader <- c("",
"", "",
paste("", kmlname, "", sep = ""), paste("", sep = ""),
paste("",0,"",sep="") )
kmlFooter <- c("", "")
kml <- append(kml, "")
kml <- append(kml, paste("", name, "", sep = ""))
kml <- append(kml, paste("", sep = ""))
kml <- append(kml, paste("#", styleNum, "",
sep = ""))
kml <- append(kml, paste("", as.integer(visibility),
"", sep = ""))
kml <- append(kml, "")
holeFlag <- FALSE
for (i in 1:length(obj@Polygons)) {
if (!holeFlag) kml <- append(kml, "")
kml <- append(kml, "relativeToGround")
kml <- append(kml, "1")
kml <- append(kml, ifelse(obj@Polygons[[i]]@hole, "",
""))
kml <- append(kml, c("", ""))
kml <- append(kml, paste(
round(coordinates(obj@Polygons[[i]])[, 1], 3),
round(coordinates(obj@Polygons[[i]])[, 2], 3), 0,sep = ","))
kml <- append(kml, c("", ""))
kml <- append(kml, ifelse(obj@Polygons[[i]]@hole, "",
""))
holeFlag <- ifelse((i + 1) <= length(obj@Polygons), obj@Polygons[[i +
1]]@hole, FALSE)
if (!holeFlag) kml <- append(kml, "")
}
kml <- append(kml, "")
kml <- append(kml, "")
if (!is.null(kmlfile))
cat(paste(c(kmlHeader, kml, kmlFooter), sep = "",
collapse = "\n"), "\n", file = kmlfile, sep = "")
else list(content = kml)
} # my.kmlPolygon
col2kmlcolor <- function(col) paste(rev(sapply(col2rgb(col,
TRUE), function(x) sprintf("%02x", x))), collapse = "")
#######################################
# MAIN PROGRAM
# cycle over variables
# * create a separate folder of polygons for each level of each variable
# * assemble into a KML file for displaying the variable
# in Google Maps or Google Earth
#######################################
for (vnum in seq(vars))
{
# select the variable that will be split into categories
V = vars[[vnum]]
# set up the KML file for this variable
kmlFile <- file(paste("postmean.",V$name,".kml", sep=""), "w")
# write the KML header
cat(my.kmlPolygon(kmlname=V$pretty.name, kmldescription="Posterior Means")$header,
file=kmlFile, sep="\n")
# write styles at the top of the file -- one per level of this variable
for (L in seq(hue)) {
lwd = 1
border = 1
kmlStyle <- ""
kmlStyle <- append(kmlStyle, paste("")
cat(kmlStyle,file=kmlFile, sep="\n")
} # for L
# assign the regions to categories, based on posterior means of selected variable
plotvar = Q[ "mean" , grep(V$name, dimnames(Q)[[2]]) ]
plotcat = as.numeric( cut(plotvar, V$brks ) )
# write outline information for each polygon in each category of this variable
for (L in seq(V$legends)) {
sel = (plotcat==L)
if (any(sel)) {
cat("",file=kmlFile,sep="\n")
cat(paste("",V$legends[L],"",sep=""),file=kmlFile,sep="\n")
out <- sapply(slot(map[sel,], "polygons"),
function(x) {
num <- as.numeric(slot(x,"ID"))+1;
my.kmlPolygon(x,
name=polyname[num],
description=lab[num],
styleNum=V$style[L]
) #my.kmlPolygon
} #function
) #sapply
cat(unlist(out), file=kmlFile, sep="\n")
} # if
cat("",file=kmlFile,sep="\n")
} # for L
# close the kml and save the file
cat(my.kmlPolygon()$footer, file=kmlFile, sep="\n")
close(kmlFile)
} # for vnum
# Note: After this program has written the .kml files, they can be converted to
# the zipped (.kmz) versions that appear on the website. Simply zip the .kml
# files and then change the extensions from .zip to .kmz