How to overlap kriging spatial prediction map on a particular area of a country map in R?

Question

I have a hourly PM10 dataset for 81 observation named "seoul032823". You can download from Here. I have performed ordinary kriging on this dataset and also got spatial map for kriging prediction. I also can show the observation data points on country map. But I cant overlap the kriging spatial prediction map on country map.

What I want to do: I want to overlap my spatial prediction map on south Korea map (not whole south korea). My area of interest is latitude 37.2N to 37.7N & Longitude 126.6E to 127.2E. That means I need to crop this area from Korea map and overlap the prediction map upon this. I also need to show the original observation data points which will follow the colour of spatial map according to concentration values. For example, I want this type of map:

My R code for kriging, and showing datapoint on korea map:

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")

I have used automap package for kriging and ggplot2 for plotting Korea map.

Answer 1

I am not too familiar with spatial analysis, so there may be issues with the projection.

First, ggplot2 works better with data.frames vs spatial objects, according to this answer citing Zev Ross. Knowing this, we can extract the kriging predictions from your kriged spatial object seoul032823_OK. The rest is relatively straightforward. You probably have to fix the longitude/latitude axes labeling and make sure the dimensions are correct on the final output. (If you do that, I can edit/append the answer to include these extra steps.)

# 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")))

Edit: OP asked for points mapped to same color scale instead of fill="yellow" defined outside the aesthetics in geom_point(). Visually, this doesn't add anything since the points blend in with the kriged background, but the code is added as requested.

Edit2: If you want the plot in the original latitude and longitude coordinates, then the different layers need to be transformed into the same coordinate system. But this transformation may result in an irregular grid that will not work for geom_tile. Solution 1: stat_summary_2d to bin and average data across the irregular grid or Solution 2: plot big square points.

#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")))