2014-09-23 71 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) 

然后我想网格单元内平均的结果(z)和产生看起来像这样(.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 2014-09-23 19:26:21

+0

感谢您的支持。我已经在g20上列出了它;对于边界点,随机分配将是我的偏好。 – 2014-09-23 19:28:07

回答

2

对此,您可以使用包sp中的over(...)函数。就我所见,根本不需要包装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*对象作为覆盖并返回与该索引的向量成x每个几何的y。在这种情况下,x是SpatialPoints对象,而y是SpatialPolygons对象。因此,over(...)标识与x中的每个点相关联的y中的多边形ID(网格单元格)。剩下的只是计算手段,将手段与原始数据框合并,并对列重新命名和重新排序,以使结果看起来像您的结果。

我调整你的代码一点,因为它没有任何意义:您创建的z值的数据帧,然后将其转换为SpatialPoints对象,该对象丢弃z值...

+0

太好了,谢谢 - 这正是我一直在寻找的。 – 2014-09-23 21:48:40