¿Recorrer 16 millones de registros con ArcPy?

13

Tengo una tabla con 8 columnas y ~ 16.7 millones de registros. Necesito ejecutar un conjunto de ecuaciones if-else en las columnas. He escrito un script usando el módulo UpdateCursor, pero después de algunos millones de registros se queda sin memoria. Me preguntaba si hay una mejor manera de procesar estos 16,7 millones de registros.

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"

ACTUALIZACIÓN # 1

Ejecuté el mismo script en una computadora con 40 gb de RAM (la computadora original tenía solo 12 gb de RAM). Se completó con éxito después de ~ 16 horas. Siento que 16 horas son demasiado largas, pero nunca he trabajado con un conjunto de datos tan grande, así que no sé qué esperar. La única adición nueva a este script es arcpy.env.parallelProcessingFactor = "100%". Estoy probando dos métodos sugeridos (1) hacer 1 millón de registros en lotes y (2) usar SearchCursor y escribir resultados en csv. Informaré sobre el progreso en breve.

ACTUALIZACIÓN # 2

¡La actualización de SearchCursor y CSV funcionó de manera brillante! No tengo los tiempos de ejecución precisos, actualizaré la publicación cuando esté en la oficina mañana, pero diría que el tiempo de ejecución aproximado es de ~ 5-6 minutos, lo cual es bastante impresionante. No lo esperaba. Estoy compartiendo mi código sin pulir, cualquier comentario y mejora son bienvenidos:

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"

ACTUALIZACIÓN # 3 Actualización final. El tiempo de ejecución total para el script es de ~ 199.6 segundos / 3.2 minutos.

cptpython
fuente
1
¿Está utilizando 64 bits (ya sea en segundo plano o servidor o Pro)?
KHibma
Olvide mencionar. Estoy ejecutando 10.4 x64 en segundo plano.
cptpython
Defensor de los demonios: ¿ha intentado ejecutarlo en primer plano o desde IDLE porque al mirar su script no necesita tener ArcMap abierto?
Hornbydd
ejecútelo como un script independiente o si conoce SQL, cargue el archivo de forma en PostgreSQL y hágalo allí
ziggy
1
Entiendo que es de código abierto, pero su proceso de aprobación lleva ~ 1-2 semanas, y esto es urgente, por lo que no creo que sea factible en este caso.
cptpython

Respuestas:

4

Puede escribir el Objectid y el resultado del cálculo (cate_2) en un archivo csv. Luego, une el csv a tu archivo original, llena un campo para preservar el resultado. De esta manera, no está actualizando la tabla con el cursor DA. Podrías usar un cursor de búsqueda.

klewis
fuente
Estaba pensando lo mismo ya que hay una discusión aquí y están hablando de conjuntos de datos aún más grandes.
Hornbydd
Gracias klewis. Esto suena prometedor. Lo intentaré junto con la sugerencia de FelixIP y una discusión interesante, aunque tendré que ejecutar esto unas pocas docenas de veces.
cptpython
¡Trabajó brillantemente! He actualizado la pregunta con el último script. ¡Gracias!
cptpython
2

Disculpas, si sigo reviviendo este viejo hilo. La idea era realizar las declaraciones if-else en el ráster combinado y luego usar el nuevo campo en Búsqueda para crear un nuevo ráster. Compliqué el problema exportando los datos como una tabla e introduje un flujo de trabajo ineficiente que fue abordado por @Alex Tereshenkov. Después de darme cuenta de lo obvio, combiné los datos en 17 consultas (1 millón cada una) como lo sugirió @FelixIP. Se tomó cada lote en promedio ~ 1.5 minutos para completar y el tiempo total de ejecución fue de ~ 23.3 minutos. Este método elimina la necesidad de uniones y creo que este método cumple mejor la tarea. Aquí hay un script revisado para referencia futura:

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()
cptpython
fuente
Entonces, solo para asegurarme de que estoy entendiendo esto correctamente. En su publicación original, dijo que cuando ejecutó esto en una computadora con 40 GB de RAM, tardó ~ 16 horas en total. Pero ahora que lo dividió en 17 lotes, y tardó ~ 23 minutos en total. ¿Es eso correcto?
ianbroad
Correcto. La primera ejecución tomó ~ 16 horas con 40 GB de RAM y la segunda ejecución tomó ~ 23 minutos + otros ~ 15 minutos para realizar Lookupy exportar el ráster con categorías recién definidas.
cptpython
Solo una nota que arcpy.env.parallelProcessingFactor = "100%"no tiene efecto en su guión. No veo ninguna herramienta allí que aproveche ese entorno.
KHibma
Estás en lo correcto. Editaré el código.
cptpython
1

Podrías intentar cambiar para usar CalculateField_management . Esto evita pasar por el uso de cursores y, por el aspecto de sus opciones para el valor de categoría, puede configurar esto como cuatro subprocesos generados secuencialmente. A medida que finaliza cada subproceso, se libera su memoria antes de comenzar el siguiente. Recibes un golpe pequeño (milisegundos) generando cada subproceso.

O, si desea mantener su enfoque actual, tenga un subproceso que tome filas x a la vez. Tenga un proceso principal para manejarlo y, como antes, sigue escalando su memoria cada vez que termina. La ventaja de hacerlo de esta manera (especialmente a través de un proceso de Python independiente) es que puedes hacer un mayor uso de todos tus núcleos como subprocesos de desove en el subproceso múltiple de Python que obtienes alrededor del GIL. Esto es posible con ArcPy y un enfoque que he usado en el pasado para hacer grandes cambios de datos. Obviamente, mantenga sus fragmentos de datos bajos; de lo contrario, se quedará sin memoria más rápido.

MappaGnosis
fuente
En mi experiencia, usar arcpy.da.UpdateCursor es mucho más rápido que arcpy.CalculateField_management. Escribí un script que se ejecuta en 55,000,000 características de construcción, fue aproximadamente 5 veces más lento con la herramienta CalculateField.
offermann
El punto es configurar cuatro subprocesos y limpiar la memoria, ya que ese es el verdadero punto crítico aquí. Como describí en el segundo párrafo, podría dividir los subprocesos por filas, pero eso requiere un poco más de administración que una sola selección.
MappaGnosis
1

La lógica de manipulación de datos se puede escribir como una instrucción UPDATE SQL utilizando una expresión CASE, que puede ejecutar utilizando GDAL / OGR, por ejemplo, a través de OSGeo4W con gdal-filegdbinstalado.

Aquí está el flujo de trabajo, que utiliza en osgeo.ogrlugar de 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)

En una tabla similar con poco más de 1 millón de registros, esta consulta tomó 18 minutos. Por lo tanto, aún podría llevar de 4 a 5 horas procesar 16 millones de registros.

Mike T
fuente
Lamentablemente, el guión es parte de un guión más grande escrito usando arcpypero agradezco la respuesta. Poco a poco estoy tratando de usar GDAL más.
cptpython
1

La actualización del código en la sección # 2 en su pregunta no muestra cómo está uniendo el .csvarchivo a la tabla original en su geodatabase de archivos. Dices que tu script tardó unos 5 minutos en ejecutarse. Esto suena justo si solo ha exportado el .csvarchivo sin hacer ninguna unión. Cuando intente devolver el .csvarchivo a ArcGIS, tendrá problemas de rendimiento.

1) No puede hacer uniones directamente desde .csvla tabla de geodatabase, porque el .csvarchivo no tiene un OID (tener un campo calculado con valores únicos no ayudará, ya que aún necesitará convertir su .csvarchivo en una tabla de geodatabase). Entonces, varios minutos para la Table To Tableherramienta GP (podría usar el in_memoryespacio de trabajo para crear una tabla temporal allí, será un poco más rápido).

2) Después de cargarlo .csven una tabla de geodatabase, querrá crear un índice en el campo donde haría la unión (en su caso, el valor de origen objectiddel .csvarchivo. Esto tomaría algunos minutos en una tabla de filas de 16 mln).

3) Entonces necesitaría usar las herramientas GP Add Joino Join FieldGP. Ninguno de los dos funcionará bien en sus mesas grandes.

4) Luego, debe hacer la Calculate Fieldherramienta GP para calcular los campos recién unidos. Muchos minutos van aquí; aún más, el cálculo del campo lleva más tiempo cuando los campos que participan en el cálculo provienen de una tabla unida.

En una palabra, no obtendrás nada cercano a los 5 minutos que mencionas. Si lo logras en una hora, me impresionaría.

Para evitar tratar con el procesamiento de grandes conjuntos de datos dentro de ArcGIS, sugiero llevar sus datos fuera de ArcGIS a un pandasmarco de datos y hacer todos sus cálculos allí. Cuando haya terminado, simplemente escriba las filas del marco de datos nuevamente en una nueva tabla de geodatabase con da.InsertCursor(o podría truncar su tabla existente y escribir sus filas en la fuente).

El código completo que he escrito para comparar esto está a continuación:

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

A continuación se muestra el resultado de Debug IO (el número informado es el número de filas en una tabla utilizada) con información sobre el tiempo de ejecución para funciones individuales:

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

Insertar una fila con da.InsertCursortoma un tiempo constante, es decir, si insertar 1 fila toma, digamos, 0.1 segundos, insertar 100 filas tomará 10 segundos. Lamentablemente, el 95% + del tiempo total de ejecución se gasta leyendo la tabla de geodatabase y luego insertando las filas nuevamente en la geodatabase.

Lo mismo es aplicable para hacer un pandasmarco de datos desde un da.SearchCursorgenerador y para calcular los campos. A medida que se duplica el número de filas en la tabla de geodatabase de origen, también lo hace el tiempo de ejecución del script anterior. Por supuesto, aún necesita usar el Python de 64 bits ya que durante la ejecución, algunas estructuras de datos más grandes se manejarán en la memoria.

Alex Tereshenkov
fuente
En realidad, iba a hacer otra pregunta que hablaría sobre las limitaciones del método que usé, porque me encontré con los problemas que ha abordado anteriormente, ¡así que gracias! Lo que estoy tratando de lograr: combinar cuatro rásteres y luego realizar una instrucción if-else basada en las columnas y escribir los resultados en una nueva columna y finalmente realizar Lookupun ráster basado en los valores de la nueva columna. Mi método tenía muchos pasos innecesarios y un flujo de trabajo ineficiente, debería haberlo mencionado en mi pregunta original. Vive y aprende. Sin embargo, probaré tu guión más adelante esta semana.
cptpython