Построение 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, я получил то, что хотел, то есть среднюю площадь пруда для каждого квадрата сетки. Мне это было очень сложно, поэтому я был бы очень благодарен, если бы кто-нибудь мог предложить решение или указать мне правильное направление.

С уважением,

Эйдан


r coordinates plot r-grid
person user2358636    schedule 19.02.2014    source источник
comment
Чтобы сделать это воспроизводимым: запустите dput на KentPonds и вставьте результат в начало вашего примера (например, KentPonds <- structure(list( ... ), class = "data.frame")). Также, IMHO, ваша пространственная методология может давать довольно большую ошибку. Сетка в 1 км, вероятно, слишком хорошее разрешение для агрегирования площади пруда только на основе местоположения центроида.   -  person Noah    schedule 19.02.2014
comment
Интересно, что в Кенте, по крайней мере, при низкой плотности пруда, где может быть ошибка, большой размер пруда сильно смещен в сторону прудов меньшего размера. Однако вы можете быть правы и может потребоваться увеличить размер сетки, но я все еще остаюсь с моей проблемой неадекватных навыков R для создания сетки.   -  person user2358636    schedule 20.02.2014

Ответы (1)


arrow_upward
1
arrow_downward

Как насчет того, чтобы просто создать несколько дополнительных столбцов для указания ссылки на сетку, а затем использовать dplyr для создания агрегированных показателей. Я также показал это в виде графика. PS Я уменьшил ваш диапазон xval и yval, чтобы элементы в предоставленном наборе данных были более заметными:

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

введите описание изображения здесь

person Troy    schedule 20.02.2014
comment
Спасибо, Трой, который отлично работал со всеми моими данными и сэкономил мне много времени. - person user2358636; 20.02.2014