5

我有一些要在 Google 地图图块上绘制的 shapefile。最有效的方法是什么?一种方法可能是使用 pkg RgoogleMaps,但是,我仍然不清楚如何执行此操作。我假设使用 PlotonStaticMap 与重新格式化 shapefile 数据的某种组合

4

2 回答 2

5

从开发人员那里,它似乎工作得很好。

  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);
于 2010-01-20T23:39:15.097 回答
0

数据来自data.gov.sg。加载绘制地图所需的包

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

使用 readOGR 读取 shapefile

 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

我们观察到坐标位于不同的投影系统中。所以我们必须把它转换成网络墨卡托投影系统。我们使用 spTransform 转换 shapefile。

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    

我们注意到变换后的投影系统。让我们将它绘制在一张底图上。

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
于 2016-11-26T14:26:28.173 回答