将GIS点与ESRI形状文件中的多边形匹配

Posted 2021-04-18

我有一个ESRI形状文件，其中包含大约150个连续地理区域（多边形）的库，它们共同构成一个地理区域。我还有一个包含60,000个事件的csv文件，每个事件都有一组唯一的x，y点坐标。 csv文件中的每个事件都发生在shape文件中150个多边形中的一个（且只有一个）中，但我不知道与每个记录关联哪个多边形。因此，我需要编写一个算法，找出多边形的身份，其中每个事件发生在60,000个事件中。我需要编写的算法输出将使我能够随后计算统计数据，例如特定多边形（地理区域）内发生的特定类型事件的可能性。

（当然，形状文件不仅仅是一个文件。它是一个目录，包含8个文件扩展名的文件，包括.dbf，.prj，.qix，.sbn，.sbx，.shp，.xml和.shx 。）

我没有ArcGIS许可证。我在http://resources.arcgis.com/content/geodatabases/10.0/file-gdb-api找到了File Geodatabase API，但我不确定它是否是正确的工具集，而且我也无法找到示例代码。

任何人都可以向我展示一些代码，用于查找哪些多边形（来自形状文件）中的每一个（来自外部数据源，如csv文件）都属于哪个？

此外，我需要指定我需要能够为每个事件的csv记录添加相关多边形的标识的特定代码。因此，仅仅绘制地图上的点以显示哪些多边形包含事件是不够的。我根本不需要可视化这些数据。相反，我需要的是能够将多边形id标记到csv文件中的每个事件记录，以便我可以进行后续的非直观的数值分析。

还将赞赏有关此主题的文章，教程和其他资源的链接。我想这是人们每天都在解决的一个问题，所以如果一个人知道如何找到它，就必须有一个既定的代码库。我每天都用Java编写代码，因此首选Java解决方案。但是，如果你有一个用不同语言编写的好代码示例，我可以从另一种语言中移植一些东西。

---

*编辑：* 我根据Spacedman的建议尝试了以下R代码，并收到以下错误消息：

> myCSV <- read.csv(file="myCSVFile.csv",head=TRUE,sep=",")  
> pts = SpatialPoints(myCSV)  
> ZipCodes = readShapeSpatial("path/myshapefile.shp")  
> overlay(myCSV,ZipCodes)  
Error in function (classes, fdef, mtable)  : unable to find an inherited method for function "overlay", for signature "data.frame", "SpatialPolygonsDataFrame"  
> 

请参阅下面的其他评论。

---

第二次编辑： 我最终使用的R代码是：

myCSV <- read.csv(file="myData.csv",head=TRUE,sep=",")  
pts = SpatialPoints(myCSV)  
ZipCodes = readShapeSpatial("myPath/ZipCodes.shp")  
write.csv(ZipCodes$ZIPCODE[overlay(pts,ZipCodes)], "ZipMatches.csv", quote=FALSE, row.names=FALSE)

注意：我必须使用：

summary(ZipCodes)  

找到包含邮政编码编码的字段的名称。在我运行摘要（ZipCodes）之前，脚本只是输出每个邮政编码的索引而不是邮政编码本身。

答案

用于执行点/多边形操作的java库是JTS：

http://www.vividsolutions.com/jts/JTSHome.htm

您可能需要其他东西来读取shapefile，也许这样：

http://openmap.bbn.com/doc/api/com/bbn/openmap/layer/shape/ShapeFile.html

或者可能是对OGR和GDAL的java绑定：

http://gdal.org/java/

但是，您可能只需要使用开源GIS软件包即可：Quantum GIS是我最喜欢的，但如果你想要一个基于Java的软件包，那么就有一些 - OpenJump，gvSIG（www.osgeo.org）事情）。

在R中，如果您已将点坐标读入矩阵或数据框：

> xy
            [,1]       [,2]
 [1,]  15.224275  -0.840066
 [2,]  -1.207046   7.959644
 [3,]   9.397658  17.426323
 [4,]  28.242840 -29.581008
 [5,]  18.664603  15.361146
 [6,]   0.975846   8.534910
 [7,] -10.941825  10.438541
 [8,] -10.331097  20.260005
 [9,]   8.105570   9.595169
[10,] -14.468177  14.366980

然后，使用maptools包及其依赖项：

> require(maptools)

首先从坐标创建一个SpatialPoints对象：

> pts = SpatialPoints(xy)

然后读入你的shapefile：

> africa=readShapeSpatial("/path/to/africa.shp")

现在重叠：

> overlay(pts,africa)
 [1] 12 20 39 27 10 55 22 33 40 45

这是shapefile中行号的向量。您可以很容易地在shapefile中查找属性：

> africa$CNTRY_NAME[overlay(pts,africa)]
 [1] Congo      Ghana      Niger      Lesotho    Chad       Togo      
 [7] Guinea     Mauritania Nigeria    Senegal   
61 Levels: Algeria Angola Benin Botswana Burkina Faso Burundi ... Zimbabwe

然后，如果要将矢量写入CSV文件，

> write.csv(africa$CNTRY_NAME[overlay(pts,africa)],file="out.csv")

生产：

"","x"
"1","Congo"
"2","Ghana"
"3","Niger"
"4","Lesotho"
"5","Chad"
"6","Togo"
"7","Guinea"
"8","Mauritania"
"9","Nigeria"
"10","Senegal"

引号用引号分隔，标题为'x'。这些东西可以通过write.csv的其他选项进行调整。

如果您的任何点落在多边形之外，则返回叠加向量将为“NA”值，您可能需要对此进行测试：

> if(any(is.na(overlay(pts,africa)))){stop("splash!")}
Error: splash!

这份厚礼？

另一答案

对于没有编程的GUI解决方案，您可以查看QGIS软件。

它将加载您的多边形形状文件没有任何问题，也可以处理creation of a point layer from your CSV file。

这里60k点应该不是问题 - 我在笔记本电脑上使用过更大的数据集。

获得多边形的属性就像在两个数据集上执行spatial join一样简单。

此操作的结果将是与您的输入数据相加的点图层加上连接多边形的属性。

如果涉及太多点击或者您需要经常重复此类分析，您可以考虑使用PyQGIS编写此过程的脚本。

对于其他解决方案，请查看GIS SE站点上的Fastest way to spatially join a point CSV with a polygon Shapefile问题。