首页 文章

如何使用R中的ggplot2在投影 Map 上绘制插值数据

提问于
浏览 次
4

我想使用ggplot2在投影 Map 上绘制一些插值数据，我已经在这个问题上工作了几个星期 . 希望有人可以帮助我，非常感谢 . shapefile和数据可以在https://www.dropbox.com/s/8wfgf8207dbh79r/gpr_000b11a_e.zip?dl=0和https://www.dropbox.com/s/9czvb35vsyf3t28/Mydata.rdata?dl=0找到 .

首先，shapefile最初使用“lon-lat”投影，我需要将其转换为Albers Equal Area（aea）投影 .

library(automap)
library(ggplot2)
library(rgdal)
load("Mydata.rdata",.GlobalEnv)
canada2<-readOGR("gpr_000b11a_e.shp", layer="gpr_000b11a_e")
g <- spTransform(canada2, CRS("+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
Borders=ggplot() +geom_polygon(data=g,aes(x=long,y=lat,group=group),fill='white',color = "black")
Borders

enter image description here

我们可以看到，我们可以正确地绘制国家 . 然后我想使用Kriging方法插入数据，代码取自Smoothing out ggplot2 map .

coordinates(Mydata)<-~longitude+latitude
proj4string(Mydata)<-CRS("+proj=longlat +datum=NAD83")
sp_mydata<-spTransform(Mydata,CRS(proj4string(g)))
Krig=autoKrige(APPT~1,sp_mydata)
interp_data = as.data.frame(Krig$krige_output)
colnames(interp_data) = c("latitude","longitude","APPT_pred","APPT_var","APPT_stdev")
interp_data=interp_data[,1:3]
ggplot(data=interp_data, aes(x=longitude, y=latitude)) + geom_tile(aes(fill=APPT_pred),color=NA)

然后我们可以看到插值数据图 .
enter image description here

最后我想结合这两个数字然后我得到以下错误： Error: Don't know how to add o to a plot

ggplot(data=interp_data, aes(x=longitude, y=latitude)) + geom_tile(aes(fill=APPT_pred),color=NA)+Borders

我的问题是：如何在 Map 上绘制插值数据，然后添加网格线（经度和纬度） . 另外，我想知道如何剪切插值数据图以适应整个加拿大 Map . 谢谢您的帮助 .

r plot ggplot2 interpolation spatial

1 回答

  • 6

    挖了一点后，我想你可能想要这个：

    Krig = autoKrige(APPT~1,sp_mydata)$krige_output
Krig = Krig[!is.na(over(Krig,as(g,"SpatialPolygons"))),]  # take only the points falling in poolygons
Krig_df = as.data.frame(Krig)
names(Krig_df) = c("APPT_pred","APPT_var","APPT_stdev","longitude","latitude")
g_fort = fortify(g)
Borders = ggplot() +
  geom_raster(data=Krig_df, aes(x=longitude, y=latitude,fill=APPT_pred))+
  geom_polygon(data=g_fort,aes(x=long,y=lat,group=group),
               fill='transparent',color = "black")+
  theme_bw()
Borders

    这使：

    ![enter image description here

    唯一的问题是您在结果 Map 中仍然有"missing"插值区域（例如，在西部） . 这是因为，从 autokrige 帮助：

    new_data：包含预测位置的sp对象 . new_data可以是点集，网格或多边形 . 不得包含NA . 如果未提供此对象，则计算默认值 . 这是通过获取input_data的凸包并在该凸包中放置大约5000个网格单元来完成的

    因此，如果不提供可行的新数据作为参数，则插值区域受输入数据集的点的凸包限制（=无外推） . 这可以使用 spsample 包中的 spsample 来解决：

    library(sp)
ptsreg <- spsample(g, 4000, type = "regular")   # Define the ouput grid - 4000 points in polygons extent
Krig = autoKrige(APPT~1,sp_mydata, new_data = ptsreg)$krige_output
Krig = Krig[!is.na(over(Krig,as(g,"SpatialPolygons"))),]  # take only the points falling in poolygons
Krig_df = as.data.frame(Krig)
names(Krig_df) = c("longitude","latitude", "APPT_pred","APPT_var","APPT_stdev")
g_fort = fortify(g)
Borders = ggplot() +
  geom_raster(data=Krig_df, aes(x=longitude, y=latitude,fill=APPT_pred))+
  geom_polygon(data=g_fort,aes(x=long,y=lat,group=group),
               fill='transparent',color = "black")+
  theme_bw()
Borders

    给出：
    enter image description here

    请注意，通过增加调用 spsample 中的插值点数，可以删除仍然具有多边形边界附近的小"holes"（因为这是一个慢操作，我只要求4000，这里）

    一个更简单的快速替代方案可能是使用包 mapview

    library(mapview)
m1 <- mapview(Krig)
m2 <- mapview(g)
m2+m1

    （您可能希望使用不太详细的多边形边界shapefile，因为这很慢）

    HTH！

    回复于

相关问题