使用池塘及其池塘面积的坐标生成 1 公里的平均池塘面积网格

Question

我有英格兰肯特郡 17,306 个池塘的北向和东向坐标，以及几乎所有池塘的面积（以平方米为单位）。我正在尝试创建一个 1 公里的网格，它给出每个网格正方形的平均池塘面积，同时为没有池塘的网格正方形提供 0 值。 我一直在寻找与我类似的问题，并找到了一个使用薄板样条算法生成英国上空网格化降雨数据的问题，然后将表面绘制到该网格上并将数据写入表格 (How to produce gridded output in R and eliminate grid squares that are not over land?)。

我已经能够在我的少量数据上使用此代码来产生类似的结果。以下是我的少量数据的示例。

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

使用所有 17,306 个池塘它会导致我的计算机崩溃，但更重要的是，我希望能够以某种方式对其进行调整，而不是通过 tps 产生的预测值，我会得到我想要的，即每个网格正方形的平均池塘面积.我一直觉得这非常困难，所以如果有人能提供解决方案或指出正确的方向，我将不胜感激。

亲切的问候，

艾丹

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