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