Ответ 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([email protected],
pred = [email protected]$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 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()
. Визуально это ничего не добавляет, так как точки смешиваются с кригированным фоном, но код добавляется в соответствии с запросом.
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([email protected], pred = [email protected]$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")))