Rのカントリーマップの特定の領域でクリギング空間予測マップを重ねる方法は？

Question

私は、 "seoul032823"という81の観測のために1時間毎のPM10データセットを持っています。 Hereからダウンロードできます。私はこのデータセットで通常のクリギングを行い、クリギング予測のための空間マップも得ました。また、観測データポイントをカントリーマップに表示することもできます。しかし、クリギング空間予測マップとカントリーマップは重ならない。

私がしたいこと：私は南朝鮮地図（南韓国全体ではありません）に私の空間予測マップを重ねたいと思います。私の関心領域は、緯度37.2N〜37.7Nです。&経度126.6E〜127.2Eです。つまり、私はこの地図を韓国の地図から切り抜いて、これに予測地図を重ねる必要があります。また、濃度値に従って空間マップの色に従う元の観測データポイントを表示する必要があります。 は例えば、私はマップのこのタイプ欲しい：クリギングため

マイRコードを、韓国の地図上のデータポイントを示す：

library(sp) 
library(gstat) 
library(automap) 
library(rgdal) 
library(e1071) 
library(dplyr) 
library(lattice) 

seoul032823 <- read.csv ("seoul032823.csv") 

#plotting the pm10 data on Korea Map 
library(ggplot2) 
library(raster) 

seoul032823 <- read.csv ("seoul032823.csv") 
skorea<- getData("GADM", country= "KOR", level=1) 
plot(skorea) 

skorea<- fortify(skorea) 
ggplot()+ 
    geom_map(data= skorea, map= skorea, aes(x=long,y=lat,map_id=id,group=group), 
      fill=NA, colour="black") + 
    geom_point(data=seoul032823, aes(x=LON, y=LAT), 
      colour= "red", alpha=0.7,na.rm=T) + 
    #scale_size(range=c(2,4))+ 
    labs(title= "PM10 Concentration in Seoul Area at South Korea", 
     x="Longitude", y= "Latitude", size="PM10(microgm/m3)")+ 
    theme(title= element_text(hjust = 0.5,vjust = 1,face= c("bold"))) 

# Reprojection 
coordinates(seoul032823) <- ~LON+LAT 
proj4string(seoul032823) <- "+proj=longlat +datum=WGS84" 
seoul032823 <- spTransform(seoul032823, CRS("+proj=utm +north +zone=52 +datum=WGS84")) 

#Creating the grid for Kriging 
LON.range <- range(as.integer(seoul032823@coords[,1 ])) + c(0,1) 
LAT.range <- range(as.integer(seoul032823@coords[,2 ])) 
seoul032823.grid <- expand.grid(LON = seq(from = LON.range[1], to = LON.range[2], by = 1500), 
           LAT = seq(from = LAT.range[1], to = LAT.range[2], by = 1500)) 
plot(seoul032823.grid) 
points(seoul032823, pch= 16,col="red") 
coordinates(seoul032823.grid)<- ~LON+LAT 
gridded(seoul032823.grid)<- T 
plot(seoul032823.grid) 
points(seoul032823, pch= 16,col="red") 

# kriging spatial prediction map 
seoul032823_OK<- autoKrige(formula = PM10~1,input_data = seoul032823, new_data = seoul032823.grid) 
pts.s <- list("sp.points", seoul032823, col = "red", pch = 16) 
automapPlot(seoul032823_OK$krige_output, "var1.pred", asp = 1, 
      sp.layout = list(pts.s), main = " Kriging Prediction") 

私は韓国のマップをプロットするためクリギングとggplot2ためautomapパッケージを使用していました。

Answer 1

私は空間解析に慣れていないので、投影に問題がある可能性があります。

最初に、ggplot2は、Zev Rossというこのanswerに従って、data.framesと空間オブジェクトの方が優れています。 これを知っていると、クリグされた空間オブジェクトseoul032823_OKからクリギング予測を抽出することができます。残りは比較的簡単です。おそらく、経度/緯度軸ラベルを修正して、最終出力で寸法が正しいことを確認する必要があります。 （あなたがそれを行う場合、私は編集/これらの追加のステップを含むように答えを追加することができます。）

# Reprojection of skorea into same coordinates as sp objects 
# Not sure if this is appropriate 
coordinates(skorea) <- ~long+lat #{sp} Convert to sp object 
proj4string(skorea) <- "+proj=longlat +datum=WGS84" #{sp} set projection attributes 
#{sp} Transform to new coordinate reference system 
skorea <- spTransform(skorea, CRS("+proj=utm +north +zone=52 +datum=WGS84")) 

#Convert spatial objects into data.frames for ggplot2 
myPoints <- data.frame(seoul032823) 
myKorea <- data.frame(skorea) 
#Extract the kriging output data into a dataframe. This is the MAIN PART! 
myKrige <- data.frame(seoul032823_OK$krige_output@coords, 
         pred = seoul032823_OK$krige_output@data$var1.pred) 
head(myKrige, 3) #Preview the data 
#  LON  LAT  pred 
#1 290853 4120600 167.8167 
#2 292353 4120600 167.5182 
#3 293853 4120600 167.1047 

#OP's original plot code, adapted here to include kriging data as geom_tile 
ggplot()+ theme_minimal() + 
    geom_tile(data = myKrige, aes(x= LON, y= LAT, fill = pred)) + 
    scale_fill_gradient2(name=bquote(atop("PM10", mu*g~m^-3)), 
         high="red", mid= "plum3", low="blue", 
         space="Lab", midpoint = median(myKrige$pred)) + 
    geom_map(data= myKorea, map= myKorea, aes(x=long,y=lat,map_id=id,group=group), 
      fill=NA, colour="black") + 
    geom_point(data=myPoints, aes(x=LON, y=LAT, fill=PM10), 
      shape=21, alpha=1,na.rm=T, size=3) + 
    coord_cartesian(xlim= LON.range, ylim= LAT.range) + 
    #scale_size(range=c(2,4))+ 
    labs(title= "PM10 Concentration in Seoul Area at South Korea", 
     x="Longitude", y= "Latitude")+ 
    theme(title= element_text(hjust = 0.5,vjust = 1,face= c("bold"))) 

編集：代わりにfill="yellow"の同じカラースケールにマッピングされたポイントを求め OPが外部に定義されました美学はgeom_point()です。視覚的には、ポイントはクリグされた背景と混ざり合っているので、これは何も追加しませんが、コードは要求どおりに追加されます。

Edit2：プロットを元の緯度と経度の座標にしたい場合は、異なるレイヤーを同じ座標系に変換する必要があります。しかし、この変換は、geom_tileで動作しない不規則なグリッドになる可能性があります。 Solution 1：stat_summary_2d、不規則なグリッドのデータを平均化して平均化するか、またはSolution 2：大きな正方形の点をプロットする。

#Reproject the krige data 
myKrige1 <- myKrige 
coordinates(myKrige1) <- ~LON+LAT 
proj4string(myKrige1) <-"+proj=utm +north +zone=52 +datum=WGS84" 
myKrige_new <- spTransform(myKrige1, CRS("+proj=longlat")) 
myKrige_new <- data.frame(myKrige_new@coords, pred = myKrige_new@data$pred) 
LON.range.new <- range(myKrige_new$LON) 
LAT.range.new <- range(myKrige_new$LAT) 

#Original seoul data have correct lat/lon data 
seoul <- read.csv ("seoul032823.csv") #Reload seoul032823 data 

#Original skorea data transformed the same was as myKrige_new 
skorea1 <- getData("GADM", country= "KOR", level=1) 
#Convert SpatialPolygonsDataFrame to dataframe (deprecated. see `broom`) 
skorea1 <- fortify(skorea1) 
coordinates(skorea1) <- ~long+lat #{sp} Convert to sp object 
proj4string(skorea1) <- "+proj=longlat +datum=WGS84" #{sp} set projection attributes 1 
#{sp} Transform to new coordinate reference system 
myKorea1 <- spTransform(skorea1, CRS("+proj=longlat")) 
myKorea1 <- data.frame(myKorea1) #Convert spatial object to data.frame for ggplot 

ggplot()+ theme_minimal() + 
    #SOLUTION 1: 
    stat_summary_2d(data=myKrige_new, aes(x = LON, y = LAT, z = pred), 
        binwidth = c(0.02,0.02)) + 
    #SOLUTION 2: Uncomment the line(s) below: 
    #geom_point(data = myKrige_new, aes(x= LON, y= LAT, fill = pred), 
    #   shape=22, size=8, colour=NA) + 
    scale_fill_gradient2(name=bquote(atop("PM10", mu*g~m^-3)), 
         high="red", mid= "plum3", low="blue", 
         space="Lab", midpoint = median(myKrige_new$pred)) + 
    geom_map(data= myKorea1, map= myKorea1, aes(x=long,y=lat,map_id=id,group=group), 
      fill=NA, colour="black") + 
    geom_point(data= seoul, aes(x=LON, y=LAT, fill=PM10), 
      shape=21, alpha=1,na.rm=T, size=3) + 
    coord_cartesian(xlim= LON.range.new, ylim= LAT.range.new) + 
    #scale_size(range=c(2,4))+ 
    labs(title= "PM10 Concentration in Seoul Area at South Korea", 
     x="Longitude", y= "Latitude")+ 
    theme(title= element_text(hjust = 0.5,vjust = 1,face= c("bold")))