ggplot2を使って投影されたマップに補間データをプロットする方法

Question

ggplot2を使って投影されたマップに補間データをプロットしたいが、私は数週間この問題を解決している。誰かが私に助けてくれることを願っています。形状ファイルとデータは、https://www.dropbox.com/s/8wfgf8207dbh79r/gpr_000b11a_e.zip?dl=0とhttps://www.dropbox.com/s/9czvb35vsyf3t28/Mydata.rdata?dl=0にあります。

まず、シェイプファイルは元々 "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 

私たちが見ることができるように、我々は正しく国をプロットすることができます。次に、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) 

次に、補間データマップが表示されます。

最後に、私は、これら2つの数字を組み合わせたいと、私は次のエラーを取得する：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 

私の質問は：どのように私は（地図上で補間データをプロットして、グリッド線を追加することができます経度と緯度）。また、カナダ地図全体に合わせて補間データマップをクリップする方法を知りました。助けてくれてありがとう。与える

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 

：

唯一の問題は、あなたがまだで、「行方不明」に補間領域を持っているということです

Answer 1

は、より多くのビットを掘り後、私はあなたがこれを望むことが推測します結果の地図（例えば、西部）。 これはautokrigeヘルプなどから、という事実によるものです：あなたが引数として実現可能なNEWDATAを提供しない場合

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

したがって、補間領域は、入力データセットの点の凸包によって制限されています（=外挿なし）。 これは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と尋ねたので）

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

HTH（これは遅いので、あなたは、あまり詳細なポリゴンの境界のシェープファイルを使用する場合があります）！