2014-02-19 51 views
1

我在英格蘭肯特郡的北部和東部座標爲17306個池塘,幾乎所有池塘的面積都是平方米。我正在嘗試創建一個1公里的網格,爲每個網格方塊提供平均池塘面積,同時爲沒有池塘的網格方塊提供0個值。 我一直在尋找類似於我的問題,並找到了一個薄板樣條算法用於在英國產生網格降雨數據,然後將曲面繪製到這個網格上並將數據寫入表格(How to produce gridded output in R and eliminate grid squares that are not over land?)。使用池塘座標和池塘區域生成1公里的平均池塘面積

我已經可以在少量的數據上使用此代碼來產生類似的結果。以下是我的一小部分數據的例子。

dput(head(KentPonds, 10))  

structure(list(Eastings = c(572745.0557, 578793.9616, 573157.8562, 
573664.2026, 572735.0952, 572738.741, 572742.0182, 572281.0791, 
572267.6893, 573673.0182), Northings = c(179326.0157, 179249.0268, 
179184.2076, 179173.6464, 179148.6766, 179123.1966, 179067.6473, 
179050.8956, 178994.7816, 178996.945), PondArea_sqm = c(448L, 
85L, 52L, 183L, 318L, 511L, 276L, 330L, 772L, 203L)), .Names = c("Eastings", 
"Northings", "PondArea_sqm"), row.names = c(NA, 10L), class = "data.frame") 

#looks like this 
Eastings Northings PondArea_sqm 
572745.1 179326.0   448 
578794.0 179249.0   85 
573157.9 179184.2   52 
573664.2 179173.6   183 
572735.1 179148.7   318 
572738.7 179123.2   511 
572742.0 179067.6   276 
572281.1 179050.9   330 
572267.7 178994.8   772 
573673.0 178996.9   203 

library(fields) 
library(maptools) 
library(gstat) 


names(KentPonds) <- c("Eastings", "Northings", "PondArea_sqm") 

fit <- Tps(cbind(KentPonds$Eastings,KentPonds$Northings),KentPonds$PondArea_sqm) 
surface(fit) 

xvals <- seq(500000, 650000, by=1000) 
yvals <- seq(115000, 190000, by=1000) 

griddf <- expand.grid(xvals, yvals) 
griddf$pred <- predict(fit, x=as.matrix(griddf)) 
write.table(griddf, file="PArea1000_Grid(1).csv", sep=",", qmethod="double") 

使用所有17306池塘它崩潰了我的電腦,但更重要的是我希望能夠以某種方式適應如此,而不是由TPS生產的預測值我會得到我想要的,即平均池塘面積是什麼爲每個網格廣場。我一直覺得這非常困難,所以如果有人能夠提供解決方案或指引我朝着正確的方向,我將非常感激。

親切的問候,

艾丹

+0

爲了使這種重複性:對'KentPonds'運行'dput',結果在貼(例如,'KentPonds < - structure(list(...),class =「data.frame」)')。另外,恕我直言,你的空間方法可能會產生相當大的錯誤。一個1公里的網格可能太精細了,只根據質心位置聚集池塘區域。 – Noah

+0

有趣的是,在肯特至少在低池塘密度,錯誤可能是高池塘大小高度偏向較小的池塘。然而,你可能是對的,可能有必要增加網格大小,但我仍然留下了我的R技能不足的問題來創建網格。 – user2358636

回答

1

怎麼樣只是創造一些額外的列來指示網格參考,然後用dplyr創建聚合措施。我也顯示它繪製。 PS我縮水的xvals和yvals範圍,從而在提供數據集中的項目會更明顯:

KentPonds<-read.table(header=T,text="Eastings Northings PondArea_sqm 
572745.1 179326.0   448 
578794.0 179249.0   85 
573157.9 179184.2   52 
573664.2 179173.6   183 
572735.1 179148.7   318 
572738.7 179123.2   511 
572742.0 179067.6   276 
572281.1 179050.9   330 
572267.7 178994.8   772 
573673.0 178996.9   203 
") 

xvals <- seq(570000, 575000, by=1000) 
yvals <- seq(178000, 181000, by=1000) 

KentPonds$x<-ceiling(KentPonds$Eastings/1000)*1000 
KentPonds$y<-ceiling(KentPonds$Northings/1000)*1000 

require(dplyr) # for aggregation 
require(ggplot2) # for plotting 

ponds_sqkm<-group_by(KentPonds,x,y) %.% 
    summarise(mean=mean(PondArea_sqm),total=sum(PondArea_sqm)) 

plotdata<-merge(ponds_sqkm,expand.grid(x=xvals,y=yvals),all.y=T) 

ggplot(plotdata) + theme_bw() + 
    geom_tile(aes(x,y,fill=mean)) + 
    scale_fill_continuous(low="white",high="red",na.value="white") 

enter image description here

+0

感謝特洛伊使用我的所有數據完美地工作,併爲我節省了很多時間。 – user2358636