Looping through 16 million records using ArcPy?

Question

I have a table with 8 columns and ~16.7 million records. I need to run a set of if-else equations on the columns. I've written a script using UpdateCursor module, but after a few million records it runs out of memory. I was wondering if there is a better way to process these 16.7 million records.

import arcpy

arcpy.TableToTable_conversion("combine_2013", "D:/mosaic.gdb", "combo_table")

c_table = "D:/mosaic.gdb/combo_table"

fields = ['dev_agg', 'herb_agg','forest_agg','wat_agg', 'cate_2']

start_time = time.time()
print "Script Started"
with arcpy.da.UpdateCursor(c_table, fields) as cursor:
    for row in cursor:
        # row's 0,1,2,3,4 = dev, herb, forest, water, category
        #classficiation water = 1; herb = 2; dev = 3; forest = 4
        if (row[3] >= 0 and row[3] > row[2]):
            row[4] = 1
        elif (row[2] >= 0 and row[2] > row[3]):
            row[4] = 4
        elif (row[1] > 180):
            row[4] = 2
        elif (row[0] > 1):
            row[4] = 3
        cursor.updateRow(row)
end_time = time.time() - start_time
print "Script Complete - " +  str(end_time) + " seconds"

UPDATE #1

I ran the same script on a computer with 40 gb RAM (the original computer had only 12 gb RAM). It successfully completed after ~16 hours. I feel that 16 hours is too long, but I've never worked with such large dataset so I don't know what to expect. The only new addition to this script is arcpy.env.parallelProcessingFactor = "100%". I'm trying two suggested methods (1) doing 1 million records in batches and (2) using SearchCursor and writing outputs to csv. I will report on the progress shortly.

UPDATE #2

The SearchCursor and CSV update worked brilliantly! I do not have the precise run times, I'll update the post when I'm in office tomorrow but I would say the approximate run time is ~5-6 minutes which is pretty impressive. I was not expecting it. I'm sharing my unpolished code any comments and improvements are welcomed:

​import arcpy, csv, time
from arcpy import env

arcpy.env.parallelProcessingFactor = "100%"

arcpy.TableToTable_conversion("D:/mosaic.gdb/combine_2013", "D:/mosaic.gdb", "combo_table")
arcpy.AddField_management("D:/mosaic.gdb/combo_table","category","SHORT")

# Table
c_table = "D:/mosaic.gdb/combo_table"
fields = ['wat_agg', 'dev_agg', 'herb_agg','forest_agg','category', 'OBJECTID']

# CSV
c_csv = open("D:/combine.csv", "w")
c_writer = csv.writer(c_csv, delimiter= ';',lineterminator='\n')
c_writer.writerow (['OID', 'CATEGORY'])
c_reader = csv.reader(c_csv)

start_time = time.time()
with arcpy.da.SearchCursor(c_table, fields) as cursor:
    for row in cursor:
        #skip file headers
        if c_reader.line_num == 1:
            continue
        # row's 0,1,2,3,4,5 = water, dev, herb, forest, category, oid
        #classficiation water = 1; dev = 2; herb = 3; ; forest = 4
        if (row[0] >= 0 and row[0] > row[3]):
            c_writer.writerow([row[5], 1])
        elif (row[1] > 1):
            c_writer.writerow([row[5], 2])
        elif (row[2] > 180):
            c_writer.writerow([row[5], 3])
        elif (row[3] >= 0 and row[3] > row[0]):
            c_writer.writerow([row[5], 4])

c_csv.close()
end_time =  time.time() - start_time
print str(end_time) + " - Seconds"

UPDATE #3 Final update. The total run time for the script is ~ 199.6 seconds / 3.2 minutes.

Answer 1

You could write the Objectid and result of the calculation (cate_2) to a csv file. Then join the csv to your original file, populate a field, to preserve the result. This way you aren't updating the table using DA cursor. You could use a Search cursor.

Answer 2

The update of the code in the section #2 in your question doesn't show how you are joining the .csv file back to the original table in your file geodatabase. You say that your script took ~5 min to run. This sounds fair if you have only exported the .csv file without doing any joins. When you will try to bring the .csv file back to ArcGIS, you will hit the performance issues.

1) You cannot do joins directly from .csv to the geodatabase table, because the .csv file doesn't have an OID (having a field calculated with unique values won't help as you still will need to convert your .csv file into a geodatabase table). So, several minutes for Table To Table GP tool (you could use the in_memory workspace to create a temp table there, will be slightly faster).

2) After you've loaded the .csv into a geodatabase table, you would want to build an index on the field you would do the join (in your case, the source objectid vaue from the .csv file. This would take some minutes on 16mln rows table.

3) Then you would need to either use the Add Join or Join Field GP tools. Neither will perform well on your large tables.

4) Afterwards, you need to do the Calculate Field GP tool to calculate the newly joined field(s). Many minutes go here; even more, the field calculation takes more time when the fields that participate in calculation are coming from a joined table.

In a word, you won't get anything close to 5min you mention. If you will make it in an hour, I would be impressed.

To avoid dealing with processing of large datasets within ArcGIS, I suggest taking your data outside of ArcGIS into a pandas data frame and do all your calculations there. When you are done, just write the data frame rows back into a new geodatabase table with da.InsertCursor (or you could truncate your existing table and write your rows into the source one).

The complete code I've written to benchmark this is below:

import time
from functools import wraps
import arcpy
import pandas as pd

def report_time(func):
    '''Decorator reporting the execution time'''
    @wraps(func)
    def wrapper(*args, **kwargs):
        start = time.time()
        result = func(*args, **kwargs)
        end = time.time()
        print(func.__name__, round(end-start,3))
        return result
    return wrapper

#----------------------------------------------------------------------
@report_time
def make_df(in_table,limit):
    columns = [f.name for f in arcpy.ListFields(in_table) if f.name != 'OBJECTID']
    cur = arcpy.da.SearchCursor(in_table,columns,'OBJECTID < {}'.format(limit))
    rows = (row for row in cur)
    df = pd.DataFrame(rows,columns=columns)
    return df

#----------------------------------------------------------------------
@report_time
def calculate_field(df):
    df.ix[(df['DataField2'] % 2 == 0), 'Category'] = 'two'
    df.ix[(df['DataField2'] % 4 == 0), 'Category'] = 'four'
    df.ix[(df['DataField2'] % 5 == 0), 'Category'] = 'five'
    df.ix[(df['DataField2'] % 10 == 0), 'Category'] = 'ten'
    df['Category'].fillna('other', inplace=True)
    return df

#----------------------------------------------------------------------
@report_time
def save_gdb_table(df,out_table):
    rows_to_write = [tuple(r[1:]) for r in df.itertuples()]
    with arcpy.da.InsertCursor(out_table,df.columns) as ins_cur:
        for row in rows_to_write:
            ins_cur.insertRow(row)

#run for tables of various sizes
for limit in [100000,500000,1000000,5000000,15000000]:
    print '{:,}'.format(limit).center(50,'-')

    in_table = r'C:\ArcGIS\scratch.gdb\BigTraffic'
    out_table = r'C:\ArcGIS\scratch.gdb\BigTrafficUpdated'
    if arcpy.Exists(out_table):
        arcpy.TruncateTable_management(out_table)

    df = make_df(in_table,limit=limit)
    df = calculate_field(df)
    save_gdb_table(df, out_table)
    print

Below is the output from the Debug IO (the number reported is the number of rows in a table used) with info on execution time for individual functions:

---------------------100,000----------------------
('make_df', 1.141)
('calculate_field', 0.042)
('save_gdb_table', 1.788)

---------------------500,000----------------------
('make_df', 4.733)
('calculate_field', 0.197)
('save_gdb_table', 8.84)

--------------------1,000,000---------------------
('make_df', 9.315)
('calculate_field', 0.392)
('save_gdb_table', 17.605)

--------------------5,000,000---------------------
('make_df', 45.371)
('calculate_field', 1.903)
('save_gdb_table', 90.797)

--------------------15,000,000--------------------
('make_df', 136.935)
('calculate_field', 5.551)
('save_gdb_table', 275.176)

Inserting a row with da.InsertCursor takes a constant time, that is, if to insert 1 row takes, say, 0.1 second, to insert 100 rows will take 10 seconds. Sadly, 95%+ of the total execution time is spent reading geodatabase table and then inserting the rows back to the geodatabase.

The same is applicable to making a pandas data frame from a da.SearchCursor generator and to calculating the field(s). As the number of rows in your source geodatabase table doubles, so does the execution time of the script above. Of course, you still need to use the 64bit Python as during the execution, some larger data structures will be handled in memory.

Answer 3

Apologies, if I keep reviving this old thread. The idea was to perform the if-else statements on the combine raster and then use the new field in Lookup to create a new raster. I complicated the problem by exporting the data as a table and introduced inefficient workflow which was addressed by @Alex Tereshenkov. After realizing the obvious, I batched the data into 17 queries (1 million each) as it was suggested by @FelixIP. It took each batch on average ~1.5 minutes to complete and total run time was ~23.3 minutes. This method eliminates need for joins and I think this method best accomplishes the task. Here's a revised script for future reference:

import arcpy, time
from arcpy import env

def cursor():
    combine = "D:/mosaic.gdb/combine_2013"
    #arcpy.AddField_management(combine,"cat_1","SHORT")
    fields = ['wat_agg', 'dev_agg', 'herb_agg','forest_agg', 'cat_1']
    batch = ['"OBJECTID" >= 1 AND "OBJECTID" <= 1000000', '"OBJECTID" >= 1000001 AND "OBJECTID" <= 2000000', '"OBJECTID" >= 2000001 AND "OBJECTID" <= 3000000', '"OBJECTID" >= 3000001 AND "OBJECTID" <= 4000000', '"OBJECTID" >= 4000001 AND "OBJECTID" <= 5000000', '"OBJECTID" >= 5000001 AND "OBJECTID" <= 6000000', '"OBJECTID" >= 6000001 AND "OBJECTID" <= 7000000', '"OBJECTID" >= 7000001 AND "OBJECTID" <= 8000000', '"OBJECTID" >= 8000001 AND "OBJECTID" <= 9000000', '"OBJECTID" >= 9000001 AND "OBJECTID" <= 10000000', '"OBJECTID" >= 10000001 AND "OBJECTID" <= 11000000', '"OBJECTID" >= 11000001 AND "OBJECTID" <= 12000000', '"OBJECTID" >= 12000001 AND "OBJECTID" <= 13000000', '"OBJECTID" >= 13000001 AND "OBJECTID" <= 14000000', '"OBJECTID" >= 14000001 AND "OBJECTID" <= 15000000', '"OBJECTID" >= 15000001 AND "OBJECTID" <= 16000000', '"OBJECTID" >= 16000001 AND "OBJECTID" <= 16757856']
    for i in batch:
        start_time = time.time()
        with arcpy.da.UpdateCursor(combine, fields, i) as cursor:
            for row in cursor:
            # row's 0,1,2,3,4,5 = water, dev, herb, forest, category
            #classficiation water = 1; dev = 2; herb = 3; ; forest = 4
                if (row[0] >= 0 and row[0] >= row[3]):
                    row[4] = 1
                elif (row[1] > 1):
                    row[4] = 2
                elif (row[2] > 180):
                    row[4] = 3
                elif (row[3] >= 0 and row[3] > row[0]):
                    row[4] = 4
                cursor.updateRow(row)
        end_time =  time.time() - start_time
        print str(end_time) + " - Seconds"

cursor()

Answer 4

You could try changing to using CalculateField_management. This avoids looping through using cursors and, by the looks of your options for the category value, you could set this up as four subprocesses spawned sequentially. As each subprocess finishes its memory is released before starting the next one. You take a small (milliseconds) hit spawning each subprocess.

Or, if you want to keep your current approach, have a subprocess that takes x-rows at a time. Have a main process to drive it and, as before, you keep scavanging your memory each time it finishes. The bonus of doing it this way (especially through a stand-alone python process) is that you can make more use of all your cores as spawning subprocesses in python's multthreading you get around the GIL. This is possible with ArcPy and an approach I've used in the past to do massive data churns. Obviously keep your chunks of data down otherwise you'll end up running out of memory faster!

Answer 5

The data manipulation logic can be written as an UPDATE SQL statement using a CASE expression, which you can execute using GDAL/OGR, e.g. via OSGeo4W with gdal-filegdb installed.

Here is the workflow, which uses osgeo.ogr instead of arcpy:

import time
from osgeo import ogr

ds = ogr.Open('D:/mosaic.gdb', 1)
if ds is None:
    raise ValueError("You don't have a 'FileGDB' driver, or the dataset doesn't exist")
sql = '''\
UPDATE combo_table SET cate_2 = CASE
    WHEN wat_agg >= 0 AND wat_agg > forest_agg THEN 1
    WHEN dev_agg > 1 THEN 2
    WHEN herb_agg > 180 THEN 3
    WHEN forest_agg >= 0 AND forest_agg > wat_agg THEN 4
    END
'''
start_time = time.time()
ds.ExecuteSQL(sql, dialect='sqlite')
ds = None  # save, close
end_time =  time.time() - start_time
print("that took %.1f seconds" % end_time)

On a similar table with just over 1 million records, this query took 18 mins. So it might still take ~4 to 5 hours to process 16 million records.