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

Question

我想使用ggplot2在投影地图上绘制一些插值数据，我已经为此问题工作了几周。希望有人能帮助我，非常感谢。形状文件和数据可以在https://www.dropbox.com/s/8wfgf8207dbh79r/gpr_000b11a_e.zip?dl=0和https://www.dropbox.com/s/9czvb35vsyf3t28/Mydata.rdata?dl=0找到。

首先，该形状文件最初使用的是“经度-纬度”投影，我需要将其转换为等面积阿尔伯斯投影（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

我们可以看到，我们可以正确地绘制国家。然后我想使用Kriging方法对数据进行插值，代码取自平滑ggplot2地图。

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)

然后我们可以看到插值数据地图。 最后我想将这两个图形结合起来，但是出现了以下错误：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 

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

- Yang Yang

---

您遇到的一个问题是interp_data对象中的纬度和经度被颠倒了。 - lbusett

为了明确，错误是因为你不能将两个使用ggplot()创建的对象相加，因此不会出现ggplot()+geom_x()+ggplot()+geom_y()。你只能添加新的图层。 - Axeman

@Axeman 非常感谢。 - Yang Yang

Answer 1

在进一步调查后，我猜想您可能需要以下内容：

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

结果如下：

问题在于，生成的地图中仍存在“缺失”的插值区域（例如，在西部）。 这是因为根据autokrige帮助文档所述：

new_data: 包含预测位置的 sp 对象。new_data 可以是一个点集、网格或多边形。不得包含 NA。如果没有提供此对象，则会计算默认值。这是通过取输入数据的凸包并在其内放置大约 5000 个网格单元来完成的。

因此，如果您不提供可行的 newdata 参数，则插值区域将受限于输入数据集的点的凸包（即不进行外推）。 您可以使用sp包中的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

这段文字的翻译如下：

得到的结果为：

请注意，你在多边形边界附近仍然有小的“孔洞”，可以通过增加调用spsample中的插值点数来消除它们（由于这是一个缓慢的操作，在此我只请求了4000个）

一个更简单快速的替代方法是使用mapview包。

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

你可能需要使用更简单的多边形边界shapefile，因为这会很慢。

希望对你有所帮助！

- lbusett

---

谢谢，地图看起来很好，但还有一个问题。一些国家的某些地区没有被“APPT”值覆盖，比如西部地区。有什么解决办法吗？非常感谢。 - Yang Yang

从autokrige帮助文件中：一个包含预测位置的sp对象。new_data可以是点集、网格或多边形。不能包含NA值。如果未提供此对象，则会计算默认值。这是通过取输入数据的凸包并在该凸包周围放置约5000个网格单元来完成的。 - lbusett

因此，如果您未提供可行的newdata作为参数，则插值区域将受到输入数据集点的凸包的限制（即无外推）。您可以使用类似于SpatialPixels或SpatialGrid的东西创建必要的数据集。我会尽快修正答案。 - lbusett

嗨。我刚刚修改了答案以解决“缺失数据”问题。 - lbusett

抱歉再次打扰您。我想知道如果我想使用颜色渐变，例如scale_fill_gradientn(colours = heat.colors(100), space = "Lab")来显示数据的渐变，我该如何做？我已经尝试了很多次，但都失败了。感谢您的时间。 - Yang Yang