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 receive a topology error notification when intersecting a point dataset (with ~5 million points), with a polygon dataset created by buffering those points by a half mile. The goal is to create a table containing the intersection of the two datasets such that I have a list of all points within that half mile radius of each starting point. I can generate effectively identical results using either an intersection or a spatial join.

My prototype of this process works fine when I work on a small subset of each dataset. When I scale up to the full dataset, the intersect operation fails with a topology error, and the spatial join fails with an out of memory error (which is plausible given the dataset size, and memory addressing limitations of a 32bit application). Is this

Much of the time I do these operations in PostGIS (successfully and easily), but on this project I'm constrained to working in ArcGIS, with the assumption that my users will have only the ArcView level of licensing. I've also done these operations in spatiallite. I'd really rather not have to pull in OGR2OGR to move the datasets to spatiallite for the processing, but can if I must.

Machine specs: Intel Core2 Quad (Q9550), Windows 7 (64bit), 8GB of ram, plenty of hard drive space

share|improve this question
    
FWIW, I've hit the exact same errors in arcpy Python scripts with large polygon file geodatabase featureclasses and Symmetrical Difference operations. I think it's a out of memory/32-bit issue. I tried repair geometry, dice, etc. and never could get my process to run. I passed it off to someone to try in GRASS. –  Chad Cooper May 9 '12 at 20:57
    
@Chad, close, not quite a 32bit issue. See my answer below with the scoop from the folks in Redlands. –  D.E.Wright May 9 '12 at 21:06
    
Did you try it using dissolve_option NONE for the buffers? –  Kirk Kuykendall May 9 '12 at 22:06
    
Consider rethinking the approach: buffering a point produces an approximate circle. It's really a polygon of somewhere between 72 and 360 points (depending on your version of ArcGIS). In the upper limit that would be 5 million * 360 = 1.8 billion points, each requiring at least 8 bytes for the two coordinates. You don't have the RAM to cache that. Even at the lower limit of 72 points/buffer, you still exceed 2 GB. Instead, work with the original points rather than their buffers. –  whuber May 10 '12 at 15:32

3 Answers 3

First the Why! The issue you are seeing is related to ESRI's scratch workspace. In ArcGIS, ESRI now uses a fGDB to store data used in temporary tasks. This workspace is built with the default value of 4GB set as the maximum space to be used; the explanation I have gotten is this is done to prevent a large run-away process with the Unlimited option setting to crash a machine during a large geoprocess.

*The best workaround!*The way you resolve this is to do a JOIN using the ONLY_FID option, this will reduce the size of your working layer in between since your temp workspace in ArcGIS is limited to 4GB in size. I have even tried running my processes totally in a ArcSDE DB and still had this issue because the tables that the software creates for the temp/scratch workspace blow-up explosively.

I actually had a ESRI guy from Redlands come in-house to track this down with me, this is the reason for the why! It is by design, we put in a feature request to have this configurable but it died in the depths of Redlands.

share|improve this answer
    
Can you not simply override the scratch workspace with arcpy.env.scratchWorkspace? –  Chad Cooper May 9 '12 at 21:18
    
No, because this is the internal ESRI working setting, versus where you are storing your temp data. We first saw this creeping up on 9.3, then worsen in 9.3.1 and finally still happening at 10. That value just tells the system where to put those temp files, not what format to use. –  D.E.Wright May 9 '12 at 21:45

Here are a few suggestions that may help:

1) Perform a "Repair Geometry" on both of the datasets before attempting the intersect, or better yet, build a quick topology and correct any polygon overlay errors as needed.

2) Convert multipart features to single part features.

3) Make sure that you set a realistic/appropriate XY tolerance and XY resolution for your data using the geoprocessor.

3) Eliminate any unnecessary fields either through layers or by choosing the "ONLY_FID" option in the Join Attributes parameter (you can get the attributes later by joining via FID/OID to the original tables).

4) Play around with different formats - despite ESRI's push to try to get everyone to use the FGDB format, the lowly Shapefile format is often much faster and less error prone when performing complex spatial operations.

share|improve this answer
    
Do you know of any benchmarks to show that the shapefile results in faster computation vs the geodatabase? I've noticed the reverse where clipping, union, intersect can complete far quicker with data in a geodatabase compared to shapefiles. –  SaultDon May 9 '12 at 20:34
    
I also noticed features stored in a geodatabase draw/render quicker on screen compared to shapefiles. Storage location can also affect performance (networked vs local). –  SaultDon May 9 '12 at 20:38
    
SaultDon, our experience in our production services has data being much more responsive in a fGDB than in a SHP. I can't give you exact benchmarks but we will see a 15%-25% improvement depending on the data being used. As we work with datasets for the entire State of Washington, this can add up. –  D.E.Wright May 9 '12 at 21:15
    
In my personal experience, I've had certain overlay (intersection, union etc.) geoprocessing operations that produce geometry/topology errors when working with the FGDB format and function properly , without errors in SHP format. –  Brent Edwards May 10 '12 at 13:11
    
Generally speaking, I develop all of my geoprocessing scripts to write a new, local FGDB to perform all processing and then push the results to the desired format (SDE, PGDB, SHP, etc.). This is almost always the fastest and most reliable method. I only experiment with other formats as a last resort, when I have tried the options mentioned. –  Brent Edwards May 10 '12 at 13:19

If doing a JOIN using the ONLY_FID option works, and your dataset is unlikely to grow significantly, that is you best bet. If not, you might be back to the old fun of tiling your dataset, processing, and merging the results. The easiest way might be to create a grid, and for each grid cell, select the point inside it, join those with the full dataset, and append the results to the output featureclass.

This won't be particularly speedy, but should get the job done for an arbitrarily large dataset.

share|improve this answer

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.