如何在 R 中的国家地图的特定区域重叠克里金空间预测图？
How to overlap kriging spatial prediction map on a particular area of a country map in R?
我有一个名为 "seoul032823" 的 81 次观察的每小时 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")
我使用了 automap 克里金法包和 ggplot2 绘制韩国地图。
我不太熟悉空间分析，所以投影可能有问题。
首先，根据这个 answer citing Zev Ross，ggplot2 与 data.frames 相比空间对象效果更好。
知道这一点，我们可以从您的克里格空间对象 seoul032823_OK 中提取克里格预测。其余的就相对简单了。您可能必须修复 longitude/latitude 轴标签并确保最终输出的尺寸正确。 （如果你这样做，我可以 edit/append 包含这些额外步骤的答案。）
# 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")))
编辑： OP 要求映射到相同色标的点，而不是 geom_point() 中在美学之外定义的 fill="yellow"。从视觉上看，这不会添加任何内容，因为这些点与 kriged 背景融为一体，但代码是按要求添加的。
Edit2: 如果要在原经纬度坐标下出图，那么需要将不同的图层转换到同一个坐标系中。但是这种转换可能会导致不规则的网格，不适用于 geom_tile。 : stat_summary_2d to bin and average data across the irregular grid or : 绘制大方块点。
#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")))
我有一个名为 "seoul032823" 的 81 次观察的每小时 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")
我使用了 automap 克里金法包和 ggplot2 绘制韩国地图。
我不太熟悉空间分析，所以投影可能有问题。
首先，根据这个 answer citing Zev Ross，ggplot2 与 data.frames 相比空间对象效果更好。
知道这一点，我们可以从您的克里格空间对象 seoul032823_OK 中提取克里格预测。其余的就相对简单了。您可能必须修复 longitude/latitude 轴标签并确保最终输出的尺寸正确。 （如果你这样做，我可以 edit/append 包含这些额外步骤的答案。）
# 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")))
编辑： OP 要求映射到相同色标的点，而不是 geom_point() 中在美学之外定义的 fill="yellow"。从视觉上看，这不会添加任何内容，因为这些点与 kriged 背景融为一体，但代码是按要求添加的。
Edit2: 如果要在原经纬度坐标下出图，那么需要将不同的图层转换到同一个坐标系中。但是这种转换可能会导致不规则的网格，不适用于 geom_tile。 stat_summary_2d to bin and average data across the irregular grid or
#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")))