2014-09-23 17 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