आर में देश के मानचित्र के किसी विशेष क्षेत्र पर क्रिगिंग स्थानिक भविष्यवाणी मानचित्र को ओवरलैप कैसे करें?

Question

मेरे पास "seoul032823" नामक 81 अवलोकन के लिए एक घंटे का पीएम 10 डेटासेट है। आप Here से डाउनलोड कर सकते हैं। मैंने इस डेटासेट पर सामान्य क्रिगिंग किया है और क्रिगिंग भविष्यवाणी के लिए स्थानिक मानचित्र भी प्राप्त किया है। मैं देश के मानचित्र पर अवलोकन डेटा बिंदु भी दिखा सकता हूं। लेकिन मैं देश के नक्शे पर क्रिंगिंग स्थानिक भविष्यवाणी मानचित्र को ओवरलैप नहीं कर सकता।

मैं क्या करना चाहता हूं: मैं दक्षिण कोरिया मानचित्र (पूरे दक्षिण कोरिया नहीं) पर अपने स्थानिक भविष्यवाणी मानचित्र को ओवरलैप करना चाहता हूं। मेरा ब्याज क्षेत्र अक्षांश 37.2N से 37.7N & रेखांश 126.6E से 127.2E है। इसका मतलब है कि मुझे कोरिया क्षेत्र से इस क्षेत्र को फसल करने और इस पर भविष्यवाणी मानचित्र को ओवरलैप करने की आवश्यकता है। मुझे मूल अवलोकन डेटा बिंदु भी दिखाना होगा जो सांद्रता मानों के अनुसार स्थानिक मानचित्र के रंग का पालन करेंगे।

मेरे kriging के लिए आर कोड, और कोरिया के नक्शे पर डाटापॉइंट दिखा: उदाहरण के लिए, मैं नक्शा के इस प्रकार चाहते

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

मैं कोरिया मानचित्र की साजिश रचने के लिए kriging के लिए automap पैकेज और ggplot2 का इस्तेमाल किया है ।

Answer 1

मैं स्थानिक विश्लेषण से बहुत परिचित नहीं हूं, इसलिए प्रक्षेपण के साथ समस्याएं हो सकती हैं।

सबसे पहले, उद्धरण Zev Ross के अनुसार ggplot2 डेटा.फ्रेम बनाम स्थानिक वस्तुओं के साथ बेहतर काम करता है। यह जानकर, हम आपके क्रिड स्पेसियल ऑब्जेक्ट 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"))) 

संपादित करें: ओपी एक ही रंग के बजाय पैमाने fill="yellow" को मैप किया अंक के लिए कहा बाहर परिभाषित geom_point() में सौंदर्यशास्त्र। दृश्यमान रूप से, यह कुछ भी नहीं जोड़ता है क्योंकि अंक क्रोधित पृष्ठभूमि के साथ मिश्रित होते हैं, लेकिन अनुरोध के रूप में कोड जोड़ा जाता है।

संपादित 2: यदि आप मूल अक्षांश और देशांतर निर्देशांक में साजिश चाहते हैं, तो विभिन्न परतों को एक ही समन्वय प्रणाली में परिवर्तित करने की आवश्यकता है। लेकिन इस परिवर्तन के परिणामस्वरूप एक अनियमित ग्रिड हो सकता है जो geom_tile के लिए काम नहीं करेगा। Solution 1: stat_summary_2d अनियमित ग्रिड या Solution 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")))