3

I have two large feature classes (3.5 mil point and 3.8 mil poly) that I would like spatially join.

Is there a faster way to do this instead of just using the classic spatial join tool?

I'm guessing this is a repeat question just can't find it.

When joining feature classes of this size I usually just let it run for a few hours to complete.

In the future, I would like to find a way to turn these types of processes from hours to minutes.

I only pull over one unique ID field to create a foreign key and then relate the two tables.

Just looking for advice or ideas.

One to One Join Intersect

Arcpy friendly Arcgis 10.1 Processor: Intel(R) Xeon(R) CPU E5-1607 0 @ 3.00 GHz RAM: 16 GB System type: 64-bit Operating System

  • 1
    I'd use x64 background GP if you don't have it already. – Paul Nov 18 '16 at 22:10
  • There is no other way I know of unless the polygons are identical and the point is exactly in the centre. then you could programatically create the polygons from centroids and copy attributes at the same time which is faster. I know this because I built such tool for a specific purpose and it works much faster than a spatial join. (I needed to join 33 million features but all my polygons were square cells and all the points exact centres) – Jakub Sisak GeoGraphics Nov 18 '16 at 22:10
  • I am guessing OP is already using background geoprocessing since spatial join craps on that many features without it. (there is a known bug) Spatial join is slow with or without background geoprocessing. – Jakub Sisak GeoGraphics Nov 18 '16 at 22:11
  • 1
    ArcGIS Pro has Pairwise Intersect, it is supposed to be faster than regular Intersect. pro.arcgis.com/en/pro-app/tool-reference/analysis/… – klewis Nov 18 '16 at 22:49
1

Possible solution using rasters:

  • Split total extent using fishnet
  • Set environment cell size
  • a) start with first tile, i.e. n = 1
  • b) Convert polygons and points to raster
  • Calculate zonal statistics (points grid = zones, polygons grid = values)
  • Count number of records in statistics table. If less than count of points, reduce cell size and go to b), otherwise transfer data from statistics table to points table, proceed with next tile, i.e. n+=1 go to a)

I tested the script below using 2.4 million points and the same number of polygons. I used fishnet equal 5*5 tiles. It means that about 100 thousands points were processed at a time. This is a snapshot of result window:

Result window

I tried spatial join on the same datasets and when it get to 5%, I cancelled it. By this time tool was running for 8 minutes. This means that with truly massive datasets (few millions), raster approach could greatly reduce processing time.

Note that I used raster to create my points and polygons, i.e. tested polygons were the most basic squares. I expect the advantage of raster approach can be even greater with more complex shapes. I could of reduce processing time further, by using fishnet with less tiles, but didn't want to risk in_memory space corruption.

SCRIPT:

import arcpy, time, os
from arcpy import env
from arcpy.sa import ZonalStatisticsAsTable as ZS
env.overwriteoutput=True
t0=time.time()
## PARAMETERS
workGDB=r'..\Scratch.gdb'
fish_net=r'D:\Scratch\fish_net.shp'
points=r'..\Scratch.gdb\points'
polygons=r'..\Scratch.gdb\polygons'
pointsGrid=r'in_memory\pntsGR'
polygonsGrid=r'in_memory\pgonsGR'
stats=r'D:\Scratch\stats.dbf'
##allStats=workGDB+os.sep+'allStats'

G=arcpy.Geometry()
tiles=arcpy.CopyFeatures_management(fish_net, G)
nTot=0
cSize=4
aDict={}
for n,tile in enumerate(tiles):
    arcpy.AddMessage("Processing tile no %i out of %i" %(n+1,len(tiles)))
    env.extent=tile.extent
    result=arcpy.GetCount_management(points)
    nPoints=int(float(result.getOutput(0)))
    while True:
        env.cellSize = cSize
        arcpy.arcpy.PointToRaster_conversion(points, "POINTID",pointsGrid)
        arcpy.FeatureToRaster_conversion(polygons, "ID",polygonsGrid)
        outZSaT = ZS(pointsGrid, "Value",polygonsGrid,stats,"DATA", "MINIMUM")
        result=arcpy.GetCount_management(stats)
        nStats=int(float(result.getOutput(0)))
        if nStats>=nPoints:break
        cSize=float(cSize)/2
    arcpy.AddMessage("Updating lookup ....")
    arcpy.SetProgressor("step", "", 0, nPoints)
    with arcpy.da.SearchCursor(stats,("Value","MIN")) as cursor:
        for v,m in cursor:
            aDict[v]=m
            arcpy.SetProgressorPosition()
    nTot+=nPoints
    arcpy.AddMessage("Processed %i points" %nTot)

arcpy.Delete_management("in_memory")
arcpy.AddMessage("Transferring results ...")
arcpy.SetProgressor("step", "", 0, nTot)
with arcpy.da.UpdateCursor(points,("pointid","PGONID")) as cursor:
    for row in cursor:
        row[1]=aDict[row[0]]
        cursor.updateRow(row)
        arcpy.SetProgressorPosition()
arcpy.AddMessage("Total time %i sec" %int(time.time()-t0))

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.