我有一个CSV,其中包含2亿个观测值,格式如下:

id,x1,y1,x2,y2,day,color
1,"-105.4652334","39.2586939","-105.4321296","39.2236632","Monday","Black"
2,"-105.3224523","39.1323299","-105.4439944","39.3352235","Tuesday","Green"
3,"-104.4233452","39.0234355","-105.4643990","39.1223435","Wednesday","Blue"


对于每组坐标(x1 / y1和x2 / y2),我想要分配它所属的美国人口普查区或人口普查区(我在此处下载了人口普查区TIGER形状文件:ftp://ftp2.census.gov/geo/tiger/TIGER2011/TRACT/tl_2011_08_tract.zip)。因此,我需要为每个观察做两次多边形点操作。匹配必须非常准确。

最快的方法是什么,包括学习软件的时间?我可以访问具有48GB内存的计算机,以防万一这可能是一个相关的限制。

多个线程建议使用PostGIS或Spatialite(Spatialite看起来更易于使用-但它与PostGIS一样有效吗?)。如果这些是最佳选择,是否必须填充空间索引(RTree)?如果是这样,怎么做(例如使用人口普查道形状文件)?对于包含示例代码(或指向示例代码的指针)的任何建议,我将深表感谢。美国人口普查数据块数据的子样本(仅y1)(100,000点)。在我终止该过程之前,花了5个多小时。我希望可以在不到40个小时的计算时间内对整个数据集实施解决方案。

提出以前提出的问题的道歉-我已仔细阅读答案,我想知道如何实施建议。我从未使用过SQL,Python,C,并且以前只使用过一次ArcGIS –我是一个完整的初学者。

评论

40小时等于每秒将近2800个多边形点操作。在我看来,这听起来是不可能的。我不知道哪个软件(ArcGIS,PostGIS,Spatialite等)最快,但是毫无疑问需要空间索引。

如果多边形不复杂,应该没问题。从索引中获取的收益(在PostGIS中)将取决于多边形的大小。较小的多边形(较小的边界框)将对索引有所帮助。可能有可能。

1249个多边形,每个多边形约有600个点。

@Uffe Kousgaard,是的,绝对有可能。你让我尝试一下。 Se在下面回答。

勇于挑战!在某些基准测试中,SpatialLite实际上比PostGIS更快地执行,但是您必须谨慎设置RTree。我还经常发现,从“内部”运行时,ArcGIS速度较慢,而从“独立” ArcPy模块“外部”运行时,ArcGIS速度更快。

#1 楼

在我的测试中,ST_DWithin比ST_Intersects更快。这是令人惊讶的,特别是因为准备好的几何算法应该在这种情况下起作用。我认为这有可能比我在这里显示的要快得多。


我做了一些测试,发现两件事几乎使速度提高了十倍。 SATA3 ssd -disks。

,然后下面的查询花费了18秒,而不是旧笔记本电脑上的62秒。
接下来,我发现在编写不需要点表的索引之前,我完全错了。有了该索引后,ST_Intersects表现出预期的效果,事情变得非常快。
我将点表中的点数增加到100万个点,并查询:

>
在72秒内运行。
由于存在1249个多边形,因此在72秒内完成了1249000000个测试。每秒约进行17000000次测试。或每秒对所有多边形测试近14000个点。

从此测试中,您要测试的400000000点应花费大约8个小时,而不会将负载分配到多个核心上。打动我:-)



首先,要可视化结果,您可以将点几何添加到结果表中,例如在QGIS中打开它,并在import_ct字段。

其次,是的,还可以通过使用右(或左)连接来获得落在任何多边形之外的点,如下所示: br />
我做了一些测试,以验证PostGIS是否可能。

首先我不了解。每行有两个点。两个点始终都在同一多边形中吗?这样就可以对其中一个点进行计算了。如果它们可以位于两个不同的多边形中,则需要一种将一个点行连接到两个多边形的方法。

从测试中看似可行,但您可能需要一些创新的解决方案才能将负载分散到多个cpu内核上。我认为是2.2GHz,2GB RAM。如果您有48个BG RAM,我想您还有更多的CPU能力。

我所做的就是创建一个具有100000点的随机点表,如下所示:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id 
FROM imported_ct , t WHERE ST_Intersects(imported_ct.geom , t.geom);


然后添加一个如下的gid:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id 
FROM imported_ct right join t ON ST_Intersects(imported_ct.the_geom , t.geom);


然后运行: />大约需要62秒(与相同点数的ArcGIS结果进行比较)。结果是一张表格,将我的表t中的点与人口普查表中的表将其连接起来。因此,如果只需检查其中一个点就足够了,那么我的旧笔记本电脑就可以使用一个内核来完成它。 />然后您可以通过针对数据库启动多个会话并运行不同的查询来手动将负载分配给多个内核。

在我的示例中,我尝试使用50000点和两个cpu内核:
在一个数据库会话上同时运行
CREATE TABLE t AS
WITH r AS
(SELECT ST_Extent(the_geom)::geometry ext FROM imported_ct)
SELECT ST_Point(x,y) AS geom FROM 
(SELECT GENERATE_SERIES(1,100000)) s,
(SELECT ST_Xmin(ext)+(random()*(ST_Xmax(ext)-ST_Xmin(ext))) x, ST_Ymin(ext)+(random()*(ST_Ymax(ext)-ST_Ymin(ext))) y FROM r
) f;


在另一个数据库会话上同时运行
ALTER TABLE t ADD COLUMN GID SERIAL;



花费了大约36秒的时间,因此它比第一个示例要慢一些,这可能取决于同时写入光盘。但是由于bith核心在同一时间工作,所以我花的时间不超过36秒。

尝试合并表t1和t2:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE ST_Dwithin(imported_ct.the_geom , t.geom,0);


用了半秒钟。使用较新的硬件并在许多内核上分配负载,即使现实世界比测试案例要慢,这绝对应该是可能的。

可能值得注意的是该示例来自Linux(Ubuntu)。使用Windows将是另一回事。但是我还有所有其他日常应用程序正在运行,因此笔记本电脑以前的负载非常重。这样也许可以很好地模拟Windows情况,而无需打开pgadmin。

评论


我只是将.tl_2011_08_trac重命名为important_ct,因为它更易于编写。因此,只需将查询中的imported_ct更改为.tl_2011_08_trac,您就可以正常进行。

–尼克拉斯·阿文(NicklasAvén)
2012年5月2日20:41

@meer顺便说一句,不建议将template_postgis_20用作将来数据库的模板。由于您似乎拥有PostGIS 2.0,因此如果您也拥有PostgreSQL 9.1,则只需创建一个新的数据库并运行“ CREATE EXTENSION POSTGIS;”即可。

–尼克拉斯·阿文(NicklasAvén)
2012年5月2日在20:43

是的,这是我几分钟前解决的另一种错字。对于那个很抱歉。还可以尝试使用ST_Intersects版本,它应该快很多。

–尼克拉斯·阿文(NicklasAvén)
2012年5月2日在22:54

@meer并非所有点都受到影响的原因是随机点放置在矩形中,我想地图并不完全是矩形。我将在帖子中进行编辑以显示如何查看结果。

–尼克拉斯·阿文(NicklasAvén)
2012年5月3日下午6:37

@Uffe Kousgaard,是的,我想你可以这样说。它一次需要一个多边形,并通过构建边缘树来准备它。然后,它检查与准备好的多边形相对应的所有点(索引已通过重叠的bbox归类为有趣点)。

–尼克拉斯·阿文(NicklasAvén)
2012年5月3日下午6:40

#2 楼

最简单的方法可能是使用PostGIS。互联网上有一些教程,介绍如何将csv / txt点数据导入PostGIS。 Link1

我不确定PostGIS中点对点搜索的性能;它应该比ArcGIS更快。 PostGIS使用的GIST空间索引非常快。 Link2 Link3

您还可以测试MongoDB地理空间索引。但这只需要更多时间就可以开始。我相信MongoDB可能真的很快。我尚未通过多边形点搜索对其进行测试,因此无法确定。