How to make a spatial join with ArcGIS with these HUGE datasets?

  • Target: almost 3 million squares 500x500m size (the attribute table is basically empty)
  • Join feature: 12000 polylines
  • Operation: one to many
  • Match: intersect

All my attempts have failed:

  • From the Spatial Join tool in the Toolbox: after 15 min the output is a file where no join has been made (Join_Count field is 0 in all rows)
  • From the dialog box: It seems to work but my estimate is that it takes 4 days to complete - without any guarantee the output will be correct.

It works when I do a Spatial Join of a small area.

Any ideas on how to do this in a proper way?

I am using file geodatabase, ArcGIS 10.2, background processing, 32 GB RAM

(The idea is to overlap those polylines to the grid (the 500x500m squares). I will then make a dissolve to find out how many lines are in each square. Finally, I will make a density map.)

  • 1
    There must exist more efficient ways than a spatial join if your only target is to find out the amount of lines in each square? I'd do it with a Python script that clips the lines for each polygon (iterate over all polygons), count the lines that are left and write to the attribute table. Optionally create output files every few ten thousand polygons or so, just to be sure. That would also allow for batch processing more easily. – Martin Jul 16 '15 at 11:36
  • Can you elaborate on that? How can you clip the lines for each polygon? Thanks! – Manuel Frias Jul 16 '15 at 11:47
  • 1
    If you have an Advanced seat of ArcGIS Desktop, the solution to this is to Intersect the lines through the polygons, then use Frequency to compute the count of lines per polygon – Vince Jul 16 '15 at 12:04
  • 1
    If you have Spatial Analyst, you don't even need to bother with the polygons, and can just use Kernel Density to make a heatmap raster from the lines. – Vince Jul 16 '15 at 12:24
  • 2
    You should be able to achieve the Intersect/Frequency solution from @Vince using a Basic license by using Intersect/Summary Statistics (with case field) instead. – PolyGeo Jul 16 '15 at 18:46

I just took a subset of the Esri "Data & Maps" US major highways (CONUS with CLASS=1 - 11,803 features), projected it to PCS_NAD_1983_CONUS_ALBERS, and generated a 2km x 2km polygon grid in envelope {-2.6m,0.2m,+2.4m,3.2m} (2500x1500 cells = 3.75m polygons) to simulate your data.

Then I used Kernel Density to make a heatmap of 2km x 2km cells (in 0.1 minutes): heatmap

Then I intersected the lines through the polygons (10.8 minutes to generate 75,221 discrete lines), and ran Frequency to generate a count of polygon-ids (0.2 minutes to identify 55701 intersected grids). From here it was easier to populate the source feature class 'icount' with a python script (which ran for 19.7 minutes):

import arcpy

arcpy.env.workspace = r'D:\Temp\gis_se\overlay.gdb'
tablePath = 'grid_frequency_tab'
layerPath = 'grid_albers'
tableCols = ["FID_grid_albers","frequency"]

count = 0
with arcpy.da.SearchCursor(tablePath, tableCols) as cursor1:
    for row1 in cursor1:
        count = count + 1
        whereClause = """{0} = {1}""".format("Objectid_1",str(row1[0]))

        with arcpy.da.UpdateCursor(layerPath, "icount", whereClause) as cursor2:
            for row2 in cursor2:
                row2[0] = row1[1]

print str(count)+" features updated"

Then I used a query definition to select only polygons with count greater than zero, and copied the result to a new polygon feature class (0.2 min): 1:21m of US

I'd recommend you include a buffer distance in your search method (Kernel Density makes this easy), since the large number of grid cells make it difficult to find large amounts of overlap: 1:1m of NY

There are many other variants possible, including using in_memory feature class to improve query performance, and using a two DA search cursors with a dictionary to compile a polygon-id/count table in memory (then updating the polygons with a DA cursor), but the process doesn't have to take days if you use the right approach.


This is an alternative way of dealing with this problem, with a Python script that clips the lines for each polygon (iterate over all polygons), count the lines that are left and write to the attribute table.

What is presented here is pieces of a code I had lying around, and is not properly tested nor a complete script. It is configured as a script tool in Arcmap. I had to change a few variable names, I might have missed one or two.

# Reading input from user
polygons = arcpy.GetParameterAsText(0)
lines  = arcpy.GetParameterAsText(1)
outcm  = arcpy.GetParameterAsText(2)

# Create temporary FC's
tempLines = arcpy.CreateFeatureclass_management('in_memory', 'tempLines', "POLYLINE", lines)
mem_polygons = "in_memory\mem_polygons"
arcpy.CopyFeatures_management(polygons, mem_polygons)

#Clip the input file for each 50 000 polygons, and process one piece at a time

cursorFields = ('SHAPE@', 'LineCount', ...) # Whatever fields you need
with  arcpy.da.UpdateCursor(mem_polygons, cursorFields) as cmrows:
    for row in cmrows:
        arcpy.env.overwriteOutput = True
        arcpy.Clip_analysis(lines, row[0], tempLines)
        arcpy.env.overwriteOutput = False

        ## Do whatever calculations you need, eg. count lines in tempLines

        ## Write to LineCount-field
        row[1] = #your count variable 


arcpy.CopyFeatures_management(mem_polygons, outcm)

Optionally you can create output files every few ten thousand polygons or so, just to be sure. That would also allow for batch processing more easily. With 32 GB RAM you should be able to process quite a few polygons with the in_memory-space. However, you can also have this write to a regular shapefile or gdb if it becomes a memory issue.

  • 1
    This solution tackles the problem from the difficult end. You don't want to perform 3 million queries against 12 thousand lines, you want to perform 12 thousand queries against 3 million polygons. – Vince Jul 16 '15 at 12:14
  • @Vince, you are of course correct about that. This was originally designed for a somewhat different purpose. However, it shouldn't require too extensive changes to query the line layer instead of the polygon layer. – Martin Jul 16 '15 at 12:50

Similar to @Vince and @Martin, if you wanted to see how many lines intersected each polygon cell, you could:

  1. Use Feature to Line to cut the lines by the polygon boundary, then remove the lines that actually belonged to the polygon (since they're both present in the output) so you just have lines belonging to the polyline FC.
  2. Use Tabulate Intersection with the count option to summarize how many "cut" lines fell inside each polygon cell.
  • 2
    Note: Tabulate Intersection requires an Advanced (ArcInfo) license – Vince Jul 16 '15 at 16:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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