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.
1 Answers
Answers 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"))) 
0 comments:
Post a Comment