多边形中的 r 点

Posted 2023-02-23

【中文标题】多边形中的 r 点

【英文标题】：r points in polygons

【发布时间】：2010-09-27 19:30:38

【问题描述】：

我有一百万个点和一个大的形状文件——8GB——太大了，无法加载到我系统上的 R 内存中。形状文件是单层的，因此给定的x、y 最多会命中一个多边形——只要它不完全在边界上！每个多边形都标有severity - 例如1，2，3。我在具有 12GB 内存的 64 位 ubuntu 机器上使用 R。

能够将数据框“标记”到多边形severity 以便我得到带有额外列的data.frame 的最简单方法是什么，即x、y、severity？

【参考方案1】：

你只有一把锤子，并不意味着所有问题都是钉子。

将您的数据加载到 PostGIS，为您的多边形构建空间索引，并执行单个 SQL 空间叠加。将结果导出回 R。

顺便说一句，说 shapefile 是 8Gb 并不是一个非常有用的信息。形状文件至少由三个文件组成，.shp 是几何图形，.dbf 是数据库，以及连接两者的 .shx。如果您的 .dbf 为 8Gb，那么您可以通过将其替换为不同的 .dbf 来轻松读取形状本身。即使 .shp 为 8Gb，它也可能只有三个多边形，在这种情况下，简化它们可能很容易。你有多少个多边形，shapefile 的 .shp 部分有多大？

【讨论】：

很高兴您发布了 Spacedman 的答案。我只是在一些 PostGIS 文档中挖掘，以弄清楚如何做到这一点，因为我认为这可能是正确的工具。

【参考方案2】：

我认为您应该预处理数据并创建一个结构，列出网格中矩形区域的可能多边形。这样，您可以减少必须检查点的多边形，并且附加结构将适合内存，因为它只有指向多边形的指针。

这里有一张图片来说明：

您想检查黄色点在哪个多边形中。您通常会检查所有多边形，但使用网格优化（橙色线，没有绘制整个网格，只是其中一个字段）您只必须检查填充的多边形，因为它们都在网格字段内或部分在网格字段内。

类似的方法是不在内存中存储所有多边形数据，而只是多边形边界框，每个多边形只需要 2 个而不是 3 个 X/Y 对（以及指向实际多边形数据的附加指针） ，但这并没有像第一个建议那样节省空间。

【讨论】：

感谢这个 schnaader - 但你能给我一些提示，让我在 R 中这样做吗？通常对于小形状文件，我可以使用库（maptools）并将它们直接读入内存并可以访问所有内容 - 但我不知道如何管理太大而无法加载的形状文件。再次感谢。

到目前为止我还没有使用 R，所以我完全不知道如何详细说明 :) 但我认为您应该尝试自己解析文件或将其转换为您自己的文件可以自己解析，理想情况下是一些大文本文件，其中每个多边形是文件中的一行。

感谢 schnaader - 我会投赞成票，但我还没有声望！ :-)

感谢京东上传图片。

你打赌！很好的答案，图表很有帮助

【参考方案3】：

我对此很感兴趣，想知道您是否在这方面取得了任何进展。既然您提出了这个问题，我想您的计算机硬件和可用于执行此相对简单操作的软件已经有所改进，以至于解决方案（如果仍然需要！）可能非常简单，尽管可能需要很长时间才能完成处理一百万个点。您可能想尝试以下方法：

#   Load relevant libraries
library(sp)
library(maptools)
library(spatstat)

#   Read shapefile. Hopefully you have a .prj file with your .shp file
#   otherwise you need to set the proj4string argument. Don't inlcude 
#   the .shp extension in the filename. I also assume that this will
#   create a SpatialPolygonsDataFrame with the "Severity" attribute
#   attached (from your .dbf file).
myshapefile <- readShapePoly("myshapefile_without_extension",     proj4string=CRS("+proj=latlong +datum=WGS84"))


#   Read your location data in. Here I assume your data has two columns X and Y giving     locations
#   Remeber that your points need to be in the same projection as your shapefile. If they aren't
#   you should look into using spTransform() on your shapefile first.
mylocs.df <- read.table(mypoints.csv, sep=",", h=TRUE)

#   Coerce X and Y coordinates to a spatial object and set projection to be the same as
#   your shapefile (WARNING: only do this if you know your points and shapefile are in
#   the same format).
mylocs.sp <- SpatialPoints(cbind(mylocs.df$X,mylocs.df$Y),     proj4string=CRS(proj4string(myshapefile))

#   Use over() to return a dataframe equal to nrows of your mylocs.df
#   with each row corresponding to a point with the attributes from the
#   poylgon in which it fell.
severity.df <- over(mylocs.sp, myshapefile)

希望该框架可以满足您的需求。您是否可以使用您现在可用的计算机/内存来做到这一点是另一回事！

【讨论】：

嗨，Simon，谢谢您 - 内存仍然是一个问题，因为其他一些形状文件和光栅运行到大约 40gb！我有 2700 万个数据点。碰巧的是，我们通过使用 python 和 gdal 找到了一个更好的更快解决方案 - 我稍后会回答自己。

【参考方案4】：

我没有一个很好的答案，但让我抛出一个想法。你能把问题翻转过来，而不是问每个点适合哪个多边形，而是“每个多边形中有哪些点”？例如，也许您可​​以将 shapefile 分解为 2000 个县，然后逐步获取每个县并检查每个点以查看它是否在该县中。如果一个点在给定的县，那么您标记它，并在下次将其从搜索中删除。

按照同样的思路，您可以将 shapefile 分成 4 个区域。然后，您可以将单个区域加上所有点放入内存中。然后只需遍历数据 4 次。

另一个想法是使用 GIS 工具来降低 shapefile 的分辨率（节点和矢量的数量）。当然，这取决于准确性对您的用例有多重要。

【讨论】：

感谢 JD - 我会投赞成票，但我还没有声望！ :-)

【参考方案5】：

我会尝试 fastshp 包。在我的粗略测试中，它明显优于 other methods for reading shapefile。它具有专门的inside 功能，我很适合您的需求。

代码应该类似于：

shp <- read.shp(YOUR_SHP_FILE, format="polygon")) 
inside(shp, x, y)

其中 x 和 y 是坐标。

如果这不起作用，我会选择 @Spacedman 提到的 PostGIS 解决方案。

【讨论】：

这个答案+1。它非常快，但目前功能有限？另外我还没有看到 shapefile 的任何绘图方法？这是在开发吗？然而，它的速度非常快。

@SimonO101 这是一个相当新的孩子（我认为）所以不能评论未来的功能。你可以plot results using ggplot2

【参考方案6】：

回答我自己的问题......并感谢大家的帮助 - 采用的最终解决方案是使用来自 python 的 gdal，它相对容易适应光栅和形状文件。一些光栅运行到大约 40gb，一些形状文件超过 8gb - 所以它们无法在我们当时拥有的任何机器上存储在内存中（现在我可以使用 128gb 内存的机器- 但我已经搬到了新的牧场！）。 python/gdal 代码能够在 1 分钟到 40 分钟内标记 2700 万个点，具体取决于 shapefile 中多边形的大小——如果有很多小多边形，它的速度非常快——如果有一些大的（250k 点) shapefile 中的多边形，速度非常慢！但是比较一下，我们以前在 oracle 空间数据库上运行它，标记这 2700 万个点大约需要 24 小时以上，或者光栅化和标记大约需要一个小时。正如 Spacedman 建议的那样，我确实尝试在我的机器上使用带有 ssd 的 postgis，但周转时间比使用 python/gdal 慢很多，因为最终解决方案不需要将 shapefile 加载到 postgis 中。总而言之，最快的方法是使用 Python/gdal：

将形状文件和 x,y csv 复制到 ssd 修改 python 脚本的配置文件，说明文件的位置以及要标记的层 并行运行几层 - 因为它受 cpu 限制而不是 i/o 限制