如何使用R中的ggplot2在投影图上绘制插值数据 [英] How to plot interpolating data on a projected map using ggplot2 in R

问题描述

我想使用ggplot2在投影地图上绘制一些插值数据，我一直在研究这个问题几个星期。希望有人能帮助我，非常感谢。 shapefile和数据可以在

<正如我们所看到的，我们可以正确地绘制国家。然后我想使用Kriging方法插入数据，代码取自

唯一的问题是，你仍然有missin g插值区域（例如，在西部）。
这是由于从 autokrige 帮助：

new_data：包含预测位置的sp对象。 new_data可以是点集，网格或多边形。不得包含NA。如果未提供此对象，则计算默认值。这是通过获取input_data的凸包并在该凸包中放置约5000个网格单元来完成的。因此，如果您不提供可行的新数据作为参数，内插区域受到输入数据集点的凸包的限制（=无外推）。
这可以在 sp 包中使用 spsample 来解决：

 库（sp）
 ptsreg < -  spsample（g，4000，type =regular）＃定义输出网格 - 多边形范围内的4000点数
 Krig = autoKrige（APPT_1，sp_mydata，new_data = ptsreg）$ krige_output 
 Krig = Krig [！is.na（over（Krig，as（g，SpatialPolygons））），]＃take只有落在poolygons中的点
 Krig_df = as.data.frame（Krig）
名称（Krig_df）= c（经度，纬度，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（）
边界
  

给出：

请注意，您仍然有多边形边界附近的小洞可以通过增加调用 spsample （因为它是一个缓慢的操作，我只要求4000，在这里）

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

  library（mapview）
 m1 < -  mapview（Krig）
 m2 < -  mapview（g）
 m2 + m1 
  

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

HTH！

I want to plot some interpolating data on a projected map using ggplot2 and I have been working on this problem for a few weeks. Hope someone can help me, thanks a lot. The shapefile and data can be found at https://www.dropbox.com/s/8wfgf8207dbh79r/gpr_000b11a_e.zip?dl=0 and https://www.dropbox.com/s/9czvb35vsyf3t28/Mydata.rdata?dl=0 .

First, the shapefile is originally using "lon-lat" projection and I need to convert it to Albers Equal Area (aea) projection.

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

As we can see, we can plot the country correctly. Then I want to interpolate the data using Kriging method, the code is taken from 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)

Then we can see the interpolating data map.

Finally I want to combine these two figures and then I get the following error: 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 

My question is: how can I plot the interpolating data on the map and then add grid lines (longitude and latitude). Also, I wonder how to clip the interpolating data map to fit the whole Canada map. Thanks for the help.

解决方案

After digging a bit more, I guess you may want this:

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

which gives:

Only problem is that you still have "missing" interpolated areas in the resulting map (e.g., on the western part). This is due to the fact that, as from autokrige help:

new_data: A sp object containing the prediction locations. new_data can be a points set, a grid or a polygon. Must not contain NA’s. If this object is not provided a default is calculated. This is done by taking the convex hull of input_data and placing around 5000 gridcells in that convex hull

Therefore, if you do not provide a feasible newdata as argument, the interpolated area is limited by the convex hull of the points of the input dataset (= no extrapolation). This can be solved using spsample insp package:

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

which gives:

Notice that the small "holes" that you still have near polygon boundaries can be removed by increasing the number of interpolation points in the call to spsample (Since it is a slow operation I only asked for 4000, here)

A simpler quick alternative could be to use package mapview

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

(you may want to use a less detailed polygon boundaries shapefiles, since this is slow)

HTH !

这篇关于如何使用R中的ggplot2在投影图上绘制插值数据的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！