Как перекрыть карту пространственного прогнозирования кригинга на определенной области карты страны в R?

Question

У меня есть почасовой набор данных PM10 для 81 наблюдения под названием "seoul032823". Вы можете скачать здесь. Я выполнил обычный кригинг для этого набора данных, а также получил пространственную карту для прогнозирования кригинга. Я также могу показать точки данных наблюдений на карте страны. Но я не могу перекрывать карту пространственного прогнозирования кригинга на карте страны.

Что я хочу сделать: я хочу наложить карту пространственного прогнозирования на карту Южной Кореи (не всю Южную Корею). Моя область интересов - широта 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 1

Я не слишком знаком с пространственным анализом, поэтому могут быть проблемы с проекцией.

Во-первых, ggplot2 лучше работает с 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")))

Изменить: OP запросил точки, сопоставленные с той же цветовой шкалой, а не fill="yellow" определяется вне эстетики в geom_point(), Визуально это ничего не добавляет, так как точки смешиваются с фоном kriged, но код добавляется по запросу.

Edit2: если вы хотите построить график в исходных координатах широты и долготы, то разные слои необходимо преобразовать в одну и ту же систему координат. Но это преобразование может привести к неправильной сетке, которая не будет работать для geom_tile, Решение 1: stat_summary_2d складывать и усреднять данные по нерегулярной сетке или Решение 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")))