2016-04-29 151 views
0

我正在嘗試使用R來執行從伊比利亞半島收集的插值數據頻率的映射。 (像這樣的https://gis.stackexchange.com/questions/147660/strange-spatial-interpolation-results-from-ordinary-kriging無法從automap包中正確使用autokrige(R無法很好地讀取預測位置)

我的問題是,由於autokrige函數的atributte new_data中的某種錯誤,繪圖並未顯示插值數據。

https://cran.r-project.org/web/packages/automap/automap.pdf new_data: 包含預測位置的sp對象。 new_data可以是一個點集,一個網格或一個多邊形。不得包含NA。如果未提供此對象,則計算默認值。這通過獲取input_data的凸包並在該凸包中放置約5000個網格單元來完成。

我認爲問題在於R不能很好地讀取轉換爲poligons的地圖,因爲如果我避免使用new_data屬性,我會得到krigging值的propper圖。但我沒有獲得伊比利亞半島的良好形狀。 有人可以幫我嗎?我將非常感激

在這裏你可以看到我的數據:http://pastebin.com/QHjn4qjP

實際的代碼:因爲我改變我的數據座標UTM投影 現在我沒有得到錯誤消息,但最後的情節不插時,整個地圖出現一個單一的顏色:(的

setwd("C:/Users/Ruth/Dropbox/R/pruebas") 

#Libraries 

library(maps) 
library(mapdata) 
library(automap) 
library(sp) 
library(maptools) 
library(raster) 
library(rgdal) 

####################MAPA############# 

#obtain the grid of the desired country 
map<-getData('GADM', country='ESP', level=0) 

#convert the grid to projected data 
mapa.utm<-spTransform(mapa3,CRSobj =CRS(" +proj=utm +zone=29 +zone=30 +zone=31 +ellps=WGS84")) 

###############################Datos####################### 

#submit the data 
data1<-read.table("FRECUENCIASH.txt",header=T) 
head(data1) 
attach(data1) 
#convert longlat coordinates to UTM 
coordinates(data1)<-c("X","Y") 
proj4string(data1) = CRS("+proj=longlat +datum=WGS84") 
data1.utm=spTransform(data1, CRS("+proj=utm +zone=29 +zone=30 +zone=31 +ellps=WGS84 ")) 

######################Kriging interpolation ##################### 


#Performs an automatic interpolation 

kriging_result<-autoKrige(Z~1,data1.utm,mapa.utm,data_variogram = data1.utm) 

#Plot the kriging 

result1<-automapPlot(kriging_result$krige_output,"var1.pred",sp.layout = list("sp.points", data1.utm));result1 

result2<-plot(kriging_result);result2 

回答

1

你是關係到你使用未投影數據作爲輸入自動映射的事實錯誤。自動映射只能處理投影數據。谷歌搜索map projections應該給你一些背景信息。將你的數據投影到一個合適的投影,您可以使用sp包中的spTransform

事實上,它沒有newdata工作是因爲在該對象投影沒有設置,所以automap不能警告你。但是,使用latlon輸入數據的automap結果不可靠。

+0

謝謝你的回答。我試過用utm座標。現在我沒有得到這個錯誤,但是最終的圖不是插值的,整個地圖出現在一個單一的顏色中。這是我的新代碼。同比是否有這個問題的原因?非常感謝 –

+0

我的第一個猜測是,觀察結果在你想要預測的區域之外。如果觀察結果很遠,那麼克里金預測收斂於平均值,這就是你觀察到的結果。檢查座標,也許你在座標變換上做了一些錯誤,例如在轉換爲utm之前設置錯誤的起始座標系。另請注意,utm與區域一起工作,請爲您的研究區域使用正確的區域。 –

+0

我再次添加了區域,我發誓我已經寫過它們,可能是所有這些嘗試和錯誤都會抹去它們,我的錯誤。順便說一下,我使用這一行:result1 <-automapPlot(kriging_result $ krige_output,「var1.pred」,sp.layout = list(「sp.points」,data1.utm)); result1查看我的點位於何處以及它們都出現在地圖內部。所以我不知道錯誤在哪裏。並再次感謝你 –