7

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

Much of the time I do these operations in PostGIS (successfully and easily), but on this project I'm constrained to working in ArcMap, with the assumption that my users will have only the ArcView (Basic) 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

5
  • 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. Commented May 9, 2012 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
    Commented May 9, 2012 at 21:06
  • Did you try it using dissolve_option NONE for the buffers? Commented May 9, 2012 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
    Commented May 10, 2012 at 15:32
  • What is the exact error message that you are seeing?
    – PolyGeo
    Commented Dec 7, 2020 at 19:14

3 Answers 3

8

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.

4
  • Can you not simply override the scratch workspace with arcpy.env.scratchWorkspace? Commented May 9, 2012 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
    Commented May 9, 2012 at 21:45
  • I'm having an issue with a large Dissolve intermittently failing; maybe this is a workaround - make the dataset less complex. Is there a bug number we can reference? In my experience that's the best way to get things from the depths of Redlands. Commented Nov 22, 2015 at 19:13
  • @CurtisPrice Sadly no; they do not see it as a bug; more of a feature. The reason for this it to keep a machine/process from going into a runaway mode and eating all drive space. It would be possible to right a recursive model that would just exponentially increase in size.
    – D.E.Wright
    Commented Nov 23, 2015 at 16:11
2

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.

5
  • 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
    Commented May 9, 2012 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
    Commented May 9, 2012 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
    Commented May 9, 2012 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. Commented May 10, 2012 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. Commented May 10, 2012 at 13:19
0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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