【问题标题】：Split longitude and latitude columns to create grid with equal sized cells拆分经度和纬度列以创建具有相同大小的单元格的网格
【发布时间】：2017-07-31 09:11:29
【问题描述】：

我有一个具有以下坐标的数据集：

lat <- c(-5.9, 35.9, 5.13, -3.4)
lon <- c(-19.9, -6, -39.9, 9.38)

所以纬度范围是 -5.9 - 35.9，而经度范围是 -39.9 - 9.38。四个纬度/地块坐标代表区域的边界（四个角）。

我想做的是创建一个网格，将这个范围分割/分隔成 20 个（或更多）相等的单元格。我尝试了以下方法，将实例分成相等的部分并创建一个新列，为每个拆分分配一个数字 (1-20)。

df$LAT_SPLIT <- NA
df$LAT_SPLIT <- as.numeric(cut2(df$LAT, g=20))

但是，每个分割的坐标范围（度）都不相同，这会创建不同大小的网格。我的问题是，如何将上述坐标分开以创建具有相等单元格的网格，同时创建一个新列，每个单元格都分配有一个数字？

我读过不同的方法，每个单元格代表纬度变化 1 度 * 经度 30 分钟，但我不知道该怎么做。我试图更改上面的代码，以便纬度的每个度数变化都会拆分纬度列，但我也不太清楚如何做到这一点，我相信你可以使用序列？即使我可以让它工作，经度仍然会有不同的范围..

我已经在 R 中尝试过，但也欢迎任何使用 Python 的建议非常期待任何可能的解决方案，谢谢！

可重现的代码

df <- structure(list(LAT = c(35.61226, 35.34986, 35.17794, 34.60425,34.40356, 33.94433, 33.41834, 16.89566, 16.89561, 16.89561),
                     LON = c(-9.604802, -9.803048, -9.921968, -10.30782, -10.44971,-10.76656, -11.13053, -24.99776, -24.99788, -24.99773)), 
                     .Names = c("LAT","LON"), class = "data.frame", row.names = c(1L, 2L, 3L, 4L, 5L,6L, 7L, 44161L, 44162L, 44163L))

【问题讨论】：

您是在寻找 Hmisc 包中的 cut2 还是 ggplot2 中的 cut_interval？
seq(minvalue, maxvalue, length.out = 21) 将为您提供 20 个等距间隔。然后，您可以将其输入 cut(values, breaks=seq(...)) 以将您的数据分配到这些间隔中。
您想要正好 20 个单元格还是需要特定的分辨率和/或正方形大小？
@Val 感谢您的回复。我并不完全需要 20 个单元，但我希望它们都具有相同的大小。我只需要将上述坐标（网格）的范围划分为单元格，以便我可以看到每个单元格有多少坐标（纬度/经度的实例）。我的 df 由大约 45,000 个具有纬度/经度坐标的实例组成。其中一些坐标在某些区域更频繁地出现，这就是为什么我需要这些单元格，以便我可以将它们分配给不同的单元格。只是不确定最好的方法是什么。如果有人知道如何编码，将不胜感激！
@MLEN 欢迎任何可行的方法，我不太挑剔，请参阅上面的评论，它可能会使我想做的事情更清楚，谢谢

标签： python r cell geospatial latitude-longitude

【解决方案1】：

这是一种生成规则多边形网格的方法。

首先我们将 data.frame 转换为SpatialPointsDataFrame

library(sp)
dfSp <- SpatialPointsDataFrame(matrix(c(df$LON, df$LAT), nrow = nrow(df)), data = df)

之后，我们使用makegrid 创建一个规则的中心点网格，并将其转换为SpatialPointsDataFrame。

grid <- makegrid(dfSp, n = 20)
gridSp <- SpatialPointsDataFrame(grid, data = data.frame(id = rownames(grid)))

为了生成多边形，我们使用了一些raster:: 函数。首先，我们创建一个RasterLayer，然后将其转换为多边形。

library(raster)
gridSpRas <- rasterFromXYZ(gridSp)
gridPoly <- rasterToPolygons(gridSpRas, dissolve = T)

此数据还具有每个多边形的数字标识符（在本例中为 gridPoly$layer）：

str(gridPoly@data)
# 'data.frame': 30 obs. of  1 variable:
#   $ layer: num  19 20 21 22 24 14 15 16 17 18 ...

我们来看看结果：

plot(gridPoly)
points(dfSp, col = "red", pch = "+")

用例：

例如，您可以像这样计算每个多边形内的点数：

gridPoly$count <- unlist(lapply(1:length(gridPoly), 
                                function (x) {length(dfSp[gridPoly[x, ], ])}))
spplot(gridPoly, zcol = "count")

【讨论】：

@Ioki 谢谢，就像一个魅力。只要有一个额外的要求。如何将 gridPoly$count 列添加到原始 df 中，以便为属于某个单元格的每个实例分配一个数字，在本例中为 1-30。这样，我可以按“网格”列对数据进行分组以进行进一步分析。
您可以使用空间连接来完成此操作，例如用例：dfSp$count <- unlist(lapply(1:length(dfSp), function(x) {gridPoly[dfSp[x, ],]$count}))。由于您的问题是关于生成网格的，您可以将答案标记为 Q 的解决方案，让未来的用户知道这种方法是有效的。
@Ioki 将其标记为已回答，谢谢！但是，上面创建新列的代码给了我以下错误：替换有 44157 行，数据有 44155。知道为什么吗？好的，我将单元格的数量更改为 30，它给了我 42 个不同的单元格（给了我 44155 行），这很好。但是当我尝试将单元格 id 分配给坐标所属的每一行时，它会分配每行的实例数，这真的很奇怪，而且需要很长时间，我尝试将 $count 更改为 $id 但这没有用要么..
“每行的实例数”在我看来好像你有一些 length() 函数，如答案中的示例所示。出于空间连接的目的，您应该查看其他线程，例如 this 或 this。
感谢那些线程帮助了我。对于那些试图解决类似问题的人，我向多边形添加了一个 id 列，如下所示：gridPoly@data$poly.ids <- 1:nrow(gridPoly)，然后使用 spatialEco 包中的以下函数将两者结合起来：pts.poly <- point.in.poly(dfSp, gridPoly)。 @Ioki 你碰巧知道如何将单元格 ID 号添加到网格中吗？