Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

edit: maybe should mention the points feature class is local and the polygons feature class is SDE

============================================================================= 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.

Thanks for your help!

share|improve this question
1  
is your os mastermap on an oracle database? you can use the database to do the intersect which will be much faster. –  Mapperz Mar 14 at 15:10
    
Good point - or likewise with PostGIS etc. If it isn't then your data is likely to be in folders by map-tile, in which case see my answer. –  MappaGnosis Mar 14 at 15:22
    
The OSMM is in an SDE database. Although I seem some kind of arcgis desktop wizkid to my colleagues, I know virtually nothing of oracle, sde and the like. I'd love to use PostGIS but I want it finished within a week, not after a lengthy requesting/application process to get postgis installed on my machine. I might try though... –  Menno Mar 14 at 15:51
    
I highly recommend using PostGIS. If you want performance - use a parallel processing function to run portions of the query concurrently on separate CPUs. On 24 CPUs we can process 15.5M points in 14M polygons in ~30 mins. Info here (thanks to Mike Gleason): geeohspatial.blogspot.com.au/2013/12/… –  minus34 Mar 15 at 5:03

2 Answers 2

up vote 2 down vote accepted

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.

share|improve this answer
    
Thank you I'll try that. Just after posting stackexchange did come up with your question from July 2011, I saw it was very similar :) –  Menno Mar 17 at 8:49
    
@Menno I'm thinking that my "question from July 2011" that you mention, which is where I learned techniques like this, will be this one: gis.stackexchange.com/questions/12448/… –  PolyGeo Mar 17 at 9:36
    
Correct, that's the one. And thanks again, because it worked perfectly! –  Menno Mar 17 at 10:57

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.

share|improve this answer
    
This sounds totally like something I'd never thought of and I'm going to try it at once, thanks a lot! –  Menno Mar 14 at 15:52
    
Unfortunately I don't have the resources, can't use my colleagues' machines to do this. Thanks for the idea anyway! –  Menno Mar 17 at 8:47

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.