2014-02-19 49 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