# Looking for fastest solution for Point in Polygon analysis of 200 million points

Apologies for asking a question that has been asked before -- I've read through answers, and I'm left wondering how to implement the recommendations. I've never used SQL, Python, C, and have only used ArcGIS once before -- I'm a complete beginner.

I have a CSV containing 200 million observations with the following format:

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"


For each set of coordinates (x1/y1 and x2/y2), I want to assign the US Census Tract or Census Block that it falls within (I downloaded the Census tract TIGER shapefile here: ftp://ftp2.census.gov/geo/tiger/TIGER2011/TRACT/tl_2011_08_tract.zip). So, I need to do a point-in-polygon operation twice for each observation. It's important that the matches be very accurate.

What is the fastest way to do this, including time to learn the software? I have access to a computer with 48GB of Memory -- in case that might be a relevant constraint.

Several threads recommend using PostGIS or Spatialite (Spatialite looks easier to use -- but is it as efficient as PostGIS?). If those are the best options, is it imperative to populate an Spatial Index (RTree?)? If so, how does one do so (e.g. using the Census Tract Shape File)? I would be extremely grateful for any recommendations that include example code (or a pointer to example code).

My first attempt (before finding this site) consisted of using ArcGIS to do a spatial join (x1/y1 only) of subsample of the data (100,000 points) on US Census Block. That took over 5 hours before I killed the process. I'm hoping for a solution that can be implemented on the entire dataset in under 40 hours of computing time.

Thank you!

-
40 hours would equal almost 2800 point-in-polygon operations per second. It just doesn't sound possible in my mind. I have no idea which piece of software (ArcGIS, PostGIS, Spatialite etc.) is the fastest, but a spatial index is without doubt needed. – Uffe Kousgaard May 2 '12 at 5:33
Should be no problem if the polygons not to complex. The gain from the index (in PostGIS) will depend on how big the polygons are. The smaller polygons (smaller bounding boxes) the more will the indexes help. Probably it is possible. – Nicklas Avén May 2 '12 at 7:11
1249 polygons with ~ 600 points per polygon. – Uffe Kousgaard May 2 '12 at 9:25
@Uffe Kousgaard, yes it is absolutely possible. You made me give it a try. Se answer below. – Nicklas Avén May 2 '12 at 10:31
Kudos for rising to the challenge! In some bench tests SpatialLite actually performs faster than PostGIS, but you have to be careful how you set up your RTrees. I have also often found ArcGIS to be slower when running from 'inside' but faster when running with a 'stand-alone' ArcPy module 'outside'. – MappaGnosis May 2 '12 at 13:55

EDIT Oh, so many typos in one post must be some sort of record. Table names was messed up, I hope it is better now.

I also realized on the way home that something is wrong here. ST_DWithin was faster in my test than ST_Intersects. That is surprising, especially since the prepared geometry algorithm is supposed to kick in on cases like this. I think there is a chance that this will be quite a lot faster than I showed here.

EDIT 2 I did some more tests and two things almost 10-doubled the speed. First, I tried on a newer computer, but still a quite ordinary laptop, maybe except from SATA3 ssd -disks.

Then the query below took 18 seconds instead of 62 seconds on the old laptop. Next I found that I was totally wrong before when I wrote that the index on the point-table wasn't necessary. With that index in place ST_Intersects behaved as expected and things became very fast. I increased the number of points in the point-table to 1 million points and the query:

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);


runs in 72 seconds. Since there is 1249 polygons, 1249000000 tests is done in 72 seconds. That makes about 17000000 tests per second. Or testing almost 14000 points against all polygons per second.

From this test your 400000000 points to test should take about 8 hours without any trouble with distributing the load to several cores. PostGIS never stops to impress me :-)

EDIT 3 First, to visualize the result you can add the point geometry to the resulting table, open it in QGIS for instance and style it with unique values on the imported_ct field.

Second, yes, you can get also the points falling outside any polygon by using a right (or left) join like this:

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);


Stop edit

I did some tests to verify if it seems possible PostGIS.

First something I don't understand. You have two points per row. Is always both points in the same polygon? Then it is enough to do the calculations on one of the points. If they can be in two different polygons you will need a way to connect one point row to two polygons.

From the tests it seems like doable, but you might need some creative solution to spread the load over more than one cpu-core.

I tested on a 4 year old laptop with dual core centrino cpu (about 2.2GHz I think) , 2GB RAM. If you have 48 BG RAM I guess you have a lot more cpu-power too.

What I did was to create a random point table with 100000 points like this:

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;


Then running :

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);


takes about 62 seconds (compare to your ArcGIS result with the same amount of points). The result is a table connecting the points in my table t with the gid in the table with census tract.

With that speed you will do 200 mill points in about 34 hours. So, if it is enough with checking one of the point, my old laptop can do it with one core.

But if you need to check both points it might be harder.

THen you can manually distribute the load to more than one core by starting multiple sessions against the db and run different queries.

In my example with 50000 points and two cpu-cores I tried :

CREATE TABLE t1 as
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE t.gid >50000 and  ST_Dwithin(imported_ct.the_geom , t.geom,0);


on one db-session at the same time as running:

CREATE TABLE t2 as
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE t.gid <=50000 and  ST_Dwithin(imported_ct.the_geom , t.geom,0);


on another db-session.

That took about 36 seconds so it is a little bit slower than the first example probably depending on disc writing at the same time. But since bith cores are working at the same time it didn't take more than 36 seconds of my time.

To union table t1 and t2 a tried:

CREATE TABLE t3 AS
SELECT * FROM t1
UNION ALL
SELECT * FROM t2;


So, with fresher hardware and distributing the load over many cores this should absolutely be possible even if the real world will be slower than the test case.

Might be worth noting that the example is from Linux (Ubuntu). Using Windows will be another story. But I have all other daily applications running so the laptop is quite heavily loaded from before. So that might simulate the windows case quite well maybe, without opening anything but pgadmin.

HTH

Nicklas

-
I just renamed .tl_2011_08_trac to imported_ct because it was easier to write. So, just change imported_ct in my query to .tl_2011_08_trac and you should go fine. – Nicklas Avén May 2 '12 at 20:41
@meer BTW, using template_postgis_20 as anything else than a template for future databases is not recommended. Since you seems to have PostGIS 2.0, if you also have PostgreSQL 9.1 you can just create a new db and run "CREATE EXTENSION POSTGIS;" – Nicklas Avén May 2 '12 at 20:43
Yes, that was another typo that I think I fixed a few minutes ago. Sorry about that. Also try the ST_Intersects version instead, that should be quite a lot faster. – Nicklas Avén May 2 '12 at 22:54
@meer The reason not every point is affected is that the random points is placed in a rectangel and I guess the map is not exactly a rectangel. I will make an edit in the post to show how to see the result. – Nicklas Avén May 3 '12 at 6:37
@Uffe Kousgaard, Yes, I guess it you can put it that way. It takes one polygon at a time and prepares it by building a tree of the edges. Then it checks all the points (that the index has sorted out as intreresting by overlapping bboxes) against that prepared polygon. – Nicklas Avén May 3 '12 at 6:40

Probably the easiest way is with PostGIS. There are some tutorials on the internet how to import csv/txt point data into PostGIS. Link1

I am not sure about performance of point-in-polygon searches in PostGIS; it should be faster than ArcGIS. GIST spatial index that PostGIS uses is pretty fast. Link2 Link3

You could also test MongoDB geospatial index. But this requires little more time to get started. I believe that MongoDB could be really fast. I haven’t tested it with point-in-polygon searches so cannot be sure.

Hope something of this helps a bit.

-