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