1

我有英格兰肯特郡 17,306 个池塘的北向和东向坐标，以及几乎所有池塘的面积（以平方米为单位）。我正在尝试创建一个 1 公里的网格，它给出每个网格正方形的平均池塘面积，同时为没有池塘的网格正方形提供 0 值。我一直在寻找与我类似的问题，并找到了一个使用薄板样条算法生成英国上空的网格化降雨数据，然后将表面绘制到该网格上并将数据写入表格的问题（如何生成网格化输出在R中并消除不在陆地上的网格方块？）。

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

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 产生的预测值，我会得到我想要的，即每个网格正方形的平均池塘面积。我一直觉得这非常困难，所以如果有人能提供解决方案或为我指明正确的方向，我将不胜感激。

亲切的问候，

艾丹

4

1 回答 1

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

在此处输入图像描述

于 2014-02-20T09:08:43.940 回答