1`PlotPolysOnStaticMap` <-structure(function#plots polygons on map 2###This function plots/overlays polygons on a map. Typically, the polygons originate from a shapefile. 3( 4 MyMap, ##<< map image returned from e.g. \code{GetMap()} 5 polys, ##<< polygons to overlay; these can be either of class \link[PBSmapping]{PolySet} from the package PBSmapping 6### or of class \link[sp]{SpatialPolygons} from the package sp 7 col, ##<< (optional) vector of colors, one for each polygon 8 border = NULL, ##<< the color to draw the border. The default, NULL, means to use \link{par}("fg"). Use border = NA to omit borders, see \link{polygon} 9 lwd = .25, ##<< line width, see \link{par} 10 verbose = 0, ##<< level of verbosity 11 add=TRUE, ##<< start a new plot or add to an existing 12 textInPolys = NULL, ##<< text to be displayed inside polygons. 13 ... ##<< further arguments passed to \code{PlotOnStaticMap} 14){ 15 ##seealso<< \link{PlotOnStaticMap} \link{mypolygon} 16 ##details<< 17 #print(str(polys)) 18 stopifnot(class(polys)[1] == "SpatialPolygons" | class(polys)[1] == "PolySet" | class(polys)[1] == "data.frame" | class(polys)[1] == "matrix") 19 if (class(polys)[1] == "SpatialPolygons") 20 polys = SpatialToPBS(polys)$xy 21 22 Rcoords <- LatLon2XY.centered(MyMap,lat= polys[,"Y"],lon= polys[,"X"]); 23 polys.XY <- as.data.frame(polys); 24 polys.XY[,"X"] <- Rcoords$newX; 25 polys.XY[,"Y"] <- Rcoords$newY; 26 27 if ( !( "PID" %in% colnames(polys.XY)) ) 28 polys.XY[,"PID"] <- 1; 29 if ( !( "SID" %in% colnames(polys.XY)) ) 30 polys.XY[,"SID"] <- 1; 31 polys.XY[,"PIDSID"] <- apply(polys.XY[,c("PID","SID")],1,paste,collapse=":") 32 if (!add) tmp <- PlotOnStaticMap(MyMap, verbose=0, ...) 33 if (verbose>1) browser() 34 #browser() 35 if (!is.null(textInPolys)) Centers = PBSmapping::calcCentroid(polys.XY) 36if (requireNamespace("PBSmapping", quietly = TRUE) & all(c("PID","X","Y","POS") %in% colnames(polys.XY)) ) { 37 attr(polys.XY, "projection") <- NULL; 38 usr <- par('usr') 39 PBSmapping::addPolys(polys.XY,col=col, border = border, lwd = lwd, xlim =usr[1:2], ylim = usr[3:4], ...) 40 if (!is.null(textInPolys)) { 41 text(Centers[,"X"],Centers[,"Y"],textInPolys,cex=0.75, col = "blue") 42 } 43} else { 44 if (!missing(col)) { 45 polys.XY[,"col"] <- col; 46 PIDtable <- as.numeric(table(polys.XY[,"PID"])); 47 SIDtable <- as.numeric(table(polys.XY[,"PIDSID"])) 48 if (length(SIDtable)==length(col)) polys.XY[,"col"] <- rep(col, SIDtable); 49 if (length(PIDtable)==length(col)) polys.XY[,"col"] <- rep(col, PIDtable); 50 } 51 52 if ( !( "col" %in% colnames(polys.XY)) ) 53 polys.XY[,"col"] <- rgb(.1,.1,.1,.05); 54 #polys.XY[polys.XY[,"PID"] == 292,"col"] <- rgb(1,0,0,.75) 55 #tmp <- by(polys.XY[,c("X","Y","col")], polys.XY[,"PID"], mypolygon, lwd = lwd, border = border); 56 #if (0){ 57 pids = unique(polys.XY[,"PIDSID"]) 58 59 for (i in pids){ 60 jj = polys.XY[,"PIDSID"] == i; 61 xx= polys.XY[jj,]; 62 if ( ( "POS" %in% colnames(xx)) ) xx <- xx[order(xx[,"POS"]),] 63 polygon( xx[, c("X","Y")], col=xx[,"col"]); 64 if (!is.null(textInPolys)) { 65 text(Centers[i,"X"],Centers[i,"Y"],textInPolys[i], col = "blue") 66 } 67 #readLines(n=1) 68 } 69 #} 70 } 71}, ex = function(){ 72if (0){ 73 #require(PBSmapping); 74 shpFile <- paste(system.file(package = "RgoogleMaps"), "/shapes/bg11_d00.shp", sep = "") 75 #shpFile <- system.file('bg11_d00.shp', package = "RgoogleMaps"); 76 77 shp=importShapefile(shpFile,projection="LL"); 78 bb <- qbbox(lat = shp[,"Y"], lon = shp[,"X"]); 79 MyMap <- GetMap.bbox(bb$lonR, bb$latR, destfile = "DC.png"); 80 PlotPolysOnStaticMap(MyMap, shp, lwd=.5, col = rgb(0.25,0.25,0.25,0.025), add = F); 81 82 #Try an open street map: 83 84 mapOSM <- GetMap.bbox(bb$lonR, bb$latR, destfile = "DC.png", type="osm"); 85 PlotPolysOnStaticMap(mapOSM, shp, lwd=.5, col = rgb(0.75,0.25,0.25,0.15), add = F); 86 87 #North Carolina SIDS data set: 88 shpFile <- system.file("shapes/sids.shp", package="maptools"); 89 shp=importShapefile(shpFile,projection="LL"); 90 bb <- qbbox(lat = shp[,"Y"], lon = shp[,"X"]); 91 MyMap <- GetMap.bbox(bb$lonR, bb$latR, destfile = "SIDS.png"); 92 #compute regularized SID rate 93 sid <- 100*attr(shp, "PolyData")$SID74/(attr(shp, "PolyData")$BIR74+500) 94 b <- as.integer(cut(sid, quantile(sid, seq(0,1,length=8)) )); 95 b[is.na(b)] <- 1; 96 opal <- col2rgb(grey.colors(7), alpha=TRUE)/255; opal["alpha",] <- 0.2; 97 shp[,"col"] <- rgb(0.1,0.1,0.1,0.2); 98 for (i in 1:length(b)) 99 shp[shp[,"PID"] == i,"col"] <- rgb(opal[1,b[i]],opal[2,b[i]],opal[3,b[i]],opal[4,b[i]]); 100 PlotPolysOnStaticMap(MyMap, shp, lwd=.5, col = shp[,"col"], add = F); 101 102 #compare the accuracy of this plot to a Google Map overlay: 103 library(maptools); 104 qk <- SpatialPointsDataFrame(as.data.frame(shp[, c("X","Y")]), as.data.frame(shp[, c("X","Y")])) 105 sp::proj4string(qk) <- CRS("+proj=longlat"); 106 tf <- "NC.counties"; 107 SGqk <- GE_SpatialGrid(qk) 108 png(file=paste(tf, ".png", sep=""), width=SGqk$width, height=SGqk$height, 109 bg="transparent") 110 par(mar=c(0,0,0,0), xaxs="i", yaxs="i");par(mai = rep(0,4)) 111 PBSmapping::plotPolys(shp, plt=NULL) 112 dev.off() 113 maptools::kmlOverlay(SGqk, paste(tf, ".kml", sep=""), paste(tf, ".png", sep="")); 114 #This kml file can now be inspected in Google Earth or Google Maps 115 116 #or choose an aspect ratio that corresponds better to North Carolina's elongated shape: 117 MyMap <- GetMap.bbox(bb$lonR, bb$latR, destfile = "SIDS.png", size = c(640, 320), zoom = 7); 118 PlotPolysOnStaticMap(MyMap, shp, lwd=.5, col = shp[,"col"], add = F); 119 } 120}) 121 122 123