5

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

pslice
  • 503
  • 1
  • 4
  • 13

2 Answers2

5

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);
daroczig
  • 28,004
  • 7
  • 90
  • 124
pslice
  • 503
  • 1
  • 4
  • 13
  • Good solution, but you should add some context : install.packages('PBSmapping', 'RgoogleMaps', 'maptools', 'rgdal', 'ReadImages'); library(PBSmapping); library(RgoogleMaps); library(maptools); Note also that if SIDS.jpg is not working, you can try SIDS.png – Etienne Racine Jan 22 '10 at 19:37
0

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
karthikbharadwaj
  • 368
  • 1
  • 7
  • 17