How to transfer attributes to thousands of local points from ArcSDE feature class with millions of polygons?

Question

I have a set of points and a huge feature class with topographic information (OS MasterMap for who knows it). I'm talking millions of polygons huge. I want to transfer attributes of the polygons to the points where they intersect.

I've tried spatial join, but it ran out of memory. I've tried identity and it's giving me no clue at all about how long it might take, though I'm going to let it run over the weekend to see where it ends.

I then thought of making a script to take one point at a time and use that to select the underlying polygons. Doing a select layer by location with the little dialog window in ArcMap works great, it takes less than a second. However doing just that from the python window takes many minutes. I've started writing a full ArcPy script to do a select by location one by one point and then update the point's attributes with the selected polygons' attributes. Updating one point took about 15 minutes. And I've got about 8000 points to do. Unfortunately my boss won't let me sit idle and wait for a script to finish if it takes 80 days to run :(

Anyone has an idea about how this might be done more quickly, and why (why ESRI, why???) SelectByLocation works so much faster in ArcMap itself?

I should mention the points feature class is local and the polygons feature class is ArcSDE

---

SOLUTION:

Among many suggestions, the only feasible one for me, being unable to install anything or use multiple workstations at the same time, was Polygeo's solution, please see the answer below. The key thing to do is: set a processing extent to a very small area around each point (do this in a cursor), make a one-polygon copy of the large dataset using this processing extent and get the attributes from the small copy.

I'm bound to try to do a query within the oracle database itself as well, but I'll need some assistance from a colleague with that. As well I'm trying to get PostGIS installed and use that. I'll report back on these as well as soon as I have results.

Answer 1

I think you should be able to get your "Updating one point took about 15 minutes" down to a few seconds by using arcpy.env.extent.

With 8,000 points the approach I am suggesting should complete within a day (worst case) - even if you write everything to disk (which I would do in initial testing), but an in_memory workspace should trim that further.

1. arcpy.da.SearchCursor iterates your point feature class to read its coordinates and an identifier
2. Select_analysis uses the identifier to copy a one point feature class from those 8,000
3. Set arcpy.env.extent to a rectangle that is say a tenth of a metre around the coordinate
4. CopyFeatures_management copies out the OSMM layer within the Geoprocessing extent i.e. almost always a single point but if you occasionally strike a boundary then you may get a few - this should take only a second or two because I frequently use this procedure on a 3.5 million polygon cadastre
5. Intersect_analysis your one point feature class with your one (or few) polys feature class. If "transfer attributes of the polygons" is not all attributes then just reading/writing them via cursors may be used to speed this up to.
6. Append_management your 8,000 intersected point feature classes back into a single feature class or, preferably use arcpy.da.InsertCursor to do this part a lot faster.

All in all, focus on testing step 4 first - if that is taking more than a few seconds then multiplying it by 8,000 becomes an issue.

Take care to turn arcpy.env.extent back to "MAXOF" once you have finished processing.

Answer 2

If you are trying to do this using ALL the OS MasterMap data, then it is a big task indeed but also something akin to some work I have done, so it can be done.

Firstly, running ArcPy from a standalone Python script always has a hit as there is time associated with loading the libraries at the start, which is already done when you open ArcMap, but you should only take the hit once per process.

Moving on to the actual problem: Frankly I'm not surprised you ran out of memory as it would have taken a monster of a computer not to have! I would (did) subdivide such tasks and using OS map sheets or subdivision of them is a useful approach but I would certianly NOT do this one point at a time. I'd write a script that selects points in one set and polygons in another based on grid which matches the OS map sheets (100km squares may be too big memory-wise so you might have to go down to 20km or something). I would then do a spatial join as you originally tried and rinse and repeat until all the 20km tiles are completed. The trick is to get the optimal size of tile so that you have enough memory to do the task in a reasonable time with minimal disk-thrashing (especially from paging) - only you can determine that because it depends both on your data, your version of ArcGIS and the machine you're doing this on.

A further refinement, if you want this task to complete sometime this century, is to farm it out as subprocesses preferably across several machines. I found a rule of thumb to be no more subprocesses than the number of cores and if you have four or more cores, leave one free for the operating system to play with. Now, when I did this sort of work, I had access to a coder who wrote a utility to automatically farm out jobs across a network of computers for my code to crunch. I guess you don't have that luxury, so you could divide the country by sets of 20km tiles and have those as a list in a simple text file which you python code will read. Install all this on as many machines as you can get hold of with ArcGIS licences and, running subprocesses, have an environment variable which you set at run time which gives a read-offset for your tiles list so that each computer handles a unique set of tiles.