2014-09-23 3 views
0

У меня есть набор данных с точками долготы/широты и значением результата для каждого набора координат. Я хотел бы создать пространственную сетку, а затем взять среднее значение результатов для координат, которые находятся в одной и той же сетке, и создать новый фрейм данных, для которого каждой координате присваивается номер сетки и имеет усредненный результат. Например, начиная с этим кодом:Анализ данных для пространственного соединения (указывает на сетку) и создания нового набора данных в R

require(sp) 
require(raster) 

frame <- data.frame(x = c(7.5, 8.2, 8.3), y = c(1,4,4.5), z = c(10,15,30)) 

coordinates(frame) <- c("x", "y") 
proj4string(frame) <- CRS("+proj=longlat") 

grid <- GridTopology(cellcentre.offset= c(0,0), cellsize = c(2,2), cells.dim = c(5,5)) 
sg <- SpatialGrid(grid) 
poly <- as.SpatialPolygons.GridTopology(grid) 
proj4string(poly) <- CRS("+proj=longlat") 

plot(poly) 
text(coordinates(poly), labels = row.names(poly), col = "gray", cex. =.6) 
points(frame$x, frame$y, col = "blue", cex = .8) 

Я тогда хотел бы усреднить результаты (г) в пределах ячеек сетки и производит dataframe, который выглядит так (наблюдения .eg):

x y z grid grid_mean 
1 7.5 1.0 10 g20  10 
2 8.2 4.0 15 g15  22.5 
3 8.3 4.5 30 g15  22.5 

Спасибо за любую помощь.

+0

кода вы предоставили первую точку на границе между g20 и g25, а не g10. Вы уверены, что сетка - это то, что вы намереваетесь сделать? Плюс, как вы хотите иметь дело с точками на границе? – jlhoward

+0

Спасибо, что поймали это. Я перечислил его в g20; для граничных точек моим предпочтением было бы случайное присваивание. –

ответ

2

Вы можете использовать функцию over(...) в пакете sp для этого. Мне не нужен пакет raster, насколько я могу судить.

require(sp) 

frame <- data.frame(x = c(7.5, 8.2, 8.3), y = c(1,4,4.5), z = c(10,15,30)) 
points <- SpatialPoints(frame) 
proj4string(points) <- CRS("+proj=longlat") 

grid <- GridTopology(cellcentre.offset= c(0,0), cellsize = c(2,2), cells.dim = c(5,5)) 
sg <- SpatialGrid(grid) 
poly <- as.SpatialPolygons.GridTopology(grid) 
proj4string(poly) <- CRS("+proj=longlat") 

# identify grids... 
result <- data.frame(frame,grid=over(points,poly)) 
# calculate means... 
result <- merge(result,aggregate(z~grid,result,mean),by="grid") 
# rename and reorder columns to make it look like your result 
colnames(result) <- c("grid","x","y","z","grid_mean") 
result <- result[,c(2,3,4,1,5)] 
result 
#  x y z grid grid_mean 
# 1 8.2 4.0 15 15  22.5 
# 2 8.3 4.5 30 15  22.5 
# 3 7.5 1.0 10 25  10.0 

over(x,y,...) функция сравнивает два Spatial* объектов как накладки и возвращает вектор с индексом в y каждой геометрии в x. В этом случае x - объект SpatialPoints, а y - объект SpatialPolygons. Таким образом, over(...) идентифицирует идентификатор многоугольника (ячейку сетки) в y, связанный с каждой точкой в ​​x. Остальное просто вычисляет средства, объединяет средства с исходным фреймом данных и переименовывает и переупорядочивает столбцы, чтобы результат выглядел как ваш результат.

Я подправил свой код немного, потому что это не имеет смысла: вы создаете фрейм данных с Z-значениями, а затем преобразовать его в SpatialPoints объект, который отбрасывает Z-значение ...

+0

Отлично, спасибо - это именно то, что я искал. –

Смежные вопросы