我有一百万个点和一个大的形状文件——8GB——太大了，无法加载到我系统上的 R 内存中。形状文件是单层的，因此给定的x、y 最多会命中一个多边形——只要它不完全在边界上！每个多边形都标有severity - 例如1，2，3。我在具有 12GB 内存的 64 位 ubuntu 机器上使用 R。
能够将数据框“标记”到多边形severity 以便我得到带有额外列的data.frame 的最简单方法是什么，即x、y、severity？
【问题讨论】：
你只有一把锤子，并不意味着所有问题都是钉子。
将您的数据加载到 PostGIS，为您的多边形构建空间索引，并执行单个 SQL 空间叠加。将结果导出回 R。
顺便说一句，说 shapefile 是 8Gb 并不是一个非常有用的信息。形状文件至少由三个文件组成，.shp 是几何图形，.dbf 是数据库，以及连接两者的 .shx。如果您的 .dbf 为 8Gb，那么您可以通过将其替换为不同的 .dbf 来轻松读取形状本身。即使 .shp 为 8Gb，它也可能只有三个多边形，在这种情况下，简化它们可能很容易。你有多少个多边形，shapefile 的 .shp 部分有多大？
【讨论】：
我认为您应该预处理数据并创建一个结构，列出网格中矩形区域的可能多边形。这样，您可以减少必须检查点的多边形，并且附加结构将适合内存，因为它只有指向多边形的指针。
这里有一张图片来说明：
您想检查黄色点在哪个多边形中。您通常会检查所有多边形，但使用网格优化（橙色线，没有绘制整个网格，只是其中一个字段）您只必须检查填充的多边形，因为它们都在网格字段内或部分在网格字段内。
类似的方法是不在内存中存储所有多边形数据，而只是多边形边界框，每个多边形只需要 2 个而不是 3 个 X/Y 对（以及指向实际多边形数据的附加指针） ，但这并没有像第一个建议那样节省空间。
【讨论】：
我对此很感兴趣，想知道您是否在这方面取得了任何进展。既然您提出了这个问题，我想您的计算机硬件和可用于执行此相对简单操作的软件已经有所改进，以至于解决方案（如果仍然需要！）可能非常简单，尽管可能需要很长时间才能完成处理一百万个点。您可能想尝试以下方法：
# 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)
希望该框架可以满足您的需求。您是否可以使用您现在可用的计算机/内存来做到这一点是另一回事！
【讨论】：
我没有一个很好的答案，但让我抛出一个想法。你能把问题翻转过来，而不是问每个点适合哪个多边形，而是“每个多边形中有哪些点”？例如，也许您可以将 shapefile 分解为 2000 个县，然后逐步获取每个县并检查每个点以查看它是否在该县中。如果一个点在给定的县，那么您标记它，并在下次将其从搜索中删除。
按照同样的思路，您可以将 shapefile 分成 4 个区域。然后，您可以将单个区域加上所有点放入内存中。然后只需遍历数据 4 次。
另一个想法是使用 GIS 工具来降低 shapefile 的分辨率（节点和矢量的数量）。当然，这取决于准确性对您的用例有多重要。
【讨论】：
我会尝试 fastshp 包。在我的粗略测试中，它明显优于 other methods for reading shapefile。它具有专门的inside 功能，我很适合您的需求。
代码应该类似于：
shp <- read.shp(YOUR_SHP_FILE, format="polygon"))
inside(shp, x, y)
其中 x 和 y 是坐标。
如果这不起作用，我会选择 @Spacedman 提到的 PostGIS 解决方案。
【讨论】：
回答我自己的问题......并感谢大家的帮助 - 采用的最终解决方案是使用来自 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：
【讨论】：