У меня есть почасовой набор данных PM10 для 81 наблюдения с именем "seoul032823". Вы можете загрузить его с Здесь. Я выполнил обычный кригинг для этого набора данных, а также получил пространственную карту для прогнозирования кригинга. Я также могу показать точки данных наблюдений на карте страны. Но я не могу перекрыть карту пространственного прогнозирования кригинга на карте страны.
Что я хочу сделать: я хочу наложить свою карту пространственного прогнозирования на карту Южной Кореи (а не всей Южной Кореи). Область моих интересов - от 37,2 до 37,7 северной широты и от 126,6 до 127,2 восточной долготы. Это означает, что мне нужно обрезать эту область с карты Кореи и наложить на нее карту прогноза. Мне также нужно показать исходные точки данных наблюдений, которые будут соответствовать цвету пространственной карты в соответствии со значениями концентрации. Например, мне нужен такой тип карты:
Мой код 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
для построения карты Кореи.