Plotting shapefiles on top of Google map tiles [cl

问题:

I have some shapefiles I want to plot over Google Maps tiles. What's the most efficient way to do this? One path might be to use the pkg RgoogleMaps, however, it is still unclear to me how to do this. I assume using PlotonStaticMap with some combination of reformatting the shapefile data

回答1:

From the developer, it seems to work nicely.

  shpFile <- system.file("shapes/sids.shp", package="maptools");  
  shp<-importShapefile(shpFile,projection="LL");  
  bb <- qbbox(lat = shp[,"Y"], lon = shp[,"X"]);  
  MyMap <- GetMap.bbox(bb$lonR, bb$latR, destfile = "SIDS.jpg");  
  #compute regularized SID rate  
  sid <- 100*attr(shp, "PolyData")$SID74/(attr(shp, "PolyData")$BIR74+500)  
  b <- as.integer(cut(sid, quantile(sid, seq(0,1,length=8)) ));  
  b[is.na(b)] <- 1;  
  opal <- col2rgb(grey.colors(7), alpha=TRUE)/255;  opal["alpha",] <- 0.2;  
  shp[,"col"] <- rgb(0.1,0.1,0.1,0.2);  
  for (i in 1:length(b)) shp[shp[,"PID"] == i,"col"] <- rgb(opal[1,b[i]],opal[2,b[i]],opal[3,b[i]],opal[4,b[i]]);  
  PlotPolysOnStaticMap(MyMap, shp, lwd=.5, col = shp[,"col"], add = F);

回答2:

The data is from data.gov.sg. Loading the packages required to plot the maps

library(rgdal)  #for shapefiles
library(ggmap)  #plotting maps
library(ggplot2) #general plots

Reading the shapefile using readOGR

 shpfile <- readOGR(dsn = 'data/street-and-places/','StreetsandPlaces',verbose = FALSE)
head(coordinates(shpfile),5)

     coords.x1 coords.x2
[1,]  28640.40  29320.77
[2,]  29429.83  28548.69
[3,]  29224.01  28360.38
[4,]  29827.20  29664.26
[5,]  28451.00  29451.78

We observe that the coordinates are in a different projection system. So we have to convert this to web mercator projection system. We transform the shapefile using spTransform.

crsobj <- CRS("+proj=longlat +datum=WGS84")   #Web Mercator projection system
shpfile_t <- spTransform(shpfile,crsobj)      #Applying projection transformation
df <- as.data.frame(coordinates(shpfile_t))   #Converting to a data frame
head(df,5)

coords.x1<dbl>coords.x2<dbl>
1   103.8391    1.281441        
2   103.8462    1.274459        
3   103.8443    1.272756        
4   103.8497    1.284547        
5   103.8374    1.282626    

We notice the transformed projection system. Let us plot it on top of a base map.

sgmap <- get_map(location="Singapore", zoom=11,   
             maptype="roadmap", source="google") #Using Osm base map of Singapore
p <- ggmap(sgmap) +
    geom_point(data = df,
    aes(x = coords.x1, y = coords.x2),
    color = 'orange',size = 1)
p