Зацикливание 16 миллионов записей с помощью ArcPy?

13

У меня есть таблица с 8 столбцами и ~ 16,7 миллионов записей. Мне нужно запустить набор уравнений if-else для столбцов. Я написал скрипт с использованием модуля UpdateCursor, но после нескольких миллионов записей ему не хватает памяти. Мне было интересно, есть ли лучший способ обработать эти 16,7 миллиона записей.

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"

ОБНОВЛЕНИЕ № 1

Я запустил тот же сценарий на компьютере с 40 ГБ ОЗУ (исходный компьютер имел только 12 ГБ ОЗУ). Успешно завершено через ~ 16 часов. Я чувствую, что 16 часов - это слишком долго, но я никогда не работал с таким большим набором данных, поэтому не знаю, чего ожидать. Единственное новое дополнение к этому скрипту arcpy.env.parallelProcessingFactor = "100%". Я пробую два предложенных метода (1) сделать миллион записей в пакетном режиме и (2) использовать SearchCursor и записать результаты в csv. Я сообщу о прогрессе в ближайшее время.

ОБНОВЛЕНИЕ № 2

Обновление SearchCursor и CSV сработало великолепно! У меня нет точного времени выполнения, я обновлю сообщение, когда завтра нахожусь в офисе, но я бы сказал, что приблизительное время выполнения составляет ~ 5-6 минут, что довольно впечатляет. Я не ожидал этого. Я делюсь своим неполированным кодом, любые комментарии и улучшения приветствуются:

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"

ОБНОВЛЕНИЕ № 3 Окончательное обновление. Общее время выполнения сценария составляет ~ 199,6 секунд / 3,2 минуты.

cptpython
источник
1
Используете ли вы 64-битный (либо фон, либо сервер, либо Pro)?
Хибма
Забыл упомянуть. Я бегу 10,4 х64 в фоновом режиме.
cptpython
Сторонники дьяволов - пробовали ли вы запускать его на переднем плане или из IDLE, поскольку, глядя на ваш скрипт, вам не нужно открывать ArcMap?
Хорнбидд
запустите его как отдельный скрипт или, если вы знаете SQL, загрузите шейп-
файл
1
Я понимаю, что это открытый исходный код, но процесс его утверждения занимает ~ 1-2 недели, и это чувствительно ко времени, поэтому я не думаю, что это осуществимо в этом случае.
cptpython

Ответы:

4

Вы можете записать Objectid и результат расчета (cate_2) в файл CSV. Затем присоедините CSV к исходному файлу, заполните поле, чтобы сохранить результат. Таким образом, вы не обновляете таблицу с помощью курсора DA. Вы можете использовать курсор поиска.

klewis
источник
Я думал то же самое, есть обсуждение здесь , и они говорят о еще большие массивы данных.
Хорнбидд
Спасибо, Клевис. Это звучит многообещающе. Я попробую это вместе с предложением FelixIP и интересным обсуждением, хотя мне придется повторить это несколько десятков раз.
cptpython
Работал блестяще! Я обновил вопрос с последним скриптом. Благодарность!
cptpython
2

Извините, если я продолжу возрождать эту старую ветку. Идея заключалась в том, чтобы выполнить операторы if-else для объединенного растра, а затем использовать новое поле в Lookup для создания нового растра. Я усложнил проблему, экспортировав данные в виде таблицы, и представил неэффективный рабочий процесс, на который обратился @Alex Tereshenkov. Поняв очевидное, я разбил данные на 17 запросов (по 1 миллиону каждый), как это было предложено @FelixIP. В среднем на каждую партию ушло ~ 1,5 минуты, а общее время работы составило ~ 23,3 минуты. Этот метод устраняет необходимость в объединениях, и я думаю, что этот метод лучше всего решает задачу. Вот пересмотренный скрипт для дальнейшего использования:

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
источник
Итак, просто чтобы убедиться, что я правильно понимаю. В своем оригинальном сообщении вы сказали, что когда вы запускали это на компьютере с 40 ГБ ОЗУ, это заняло ~ 16 часов. Но теперь, когда вы разделили его на 17 партий, это заняло ~ 23 минуты. Это верно?
ianbroad
Верный. Первый запуск занял ~ 16 часов с 40 ГБ ОЗУ, а второй запуск занял ~ 23 минуты + еще ~ 15 минут, чтобы выполнить Lookupи экспортировать растр с новыми определенными категориями.
cptpython 13.12.16
Просто примечание, arcpy.env.parallelProcessingFactor = "100%"которое не влияет на ваш сценарий. Я не вижу никаких инструментов, которые могли бы использовать эту среду.
Хибма
Ты прав. Я буду редактировать код.
cptpython
1

Вы можете попробовать перейти на использование CalculateField_management . Это позволяет избежать циклического использования курсоров, и, судя по вашим параметрам для значения категории, вы можете настроить это как четыре подпроцесса, порожденных последовательно. По завершении каждого подпроцесса его память освобождается перед запуском следующего. Вы выполняете небольшое (миллисекунды) попадание, порождая каждый подпроцесс.

Или, если вы хотите сохранить свой текущий подход, используйте подпроцесс, который принимает по x строк одновременно. Иметь основной процесс для управления им, и, как и прежде, вы продолжаете ломать свою память каждый раз, когда она заканчивается. Преимущество такого подхода (особенно через автономный процесс Python) заключается в том, что вы можете более активно использовать все свои ядра в качестве порождающих подпроцессов в многопоточности Python, которую вы получаете вокруг GIL. Это возможно с ArcPy и подходом, который я использовал в прошлом, чтобы делать большие объемы данных. Очевидно, что ваши куски данных будут неактивными, иначе у вас закончится нехватка памяти быстрее!

MappaGnosis
источник
По моему опыту, использование arcpy.da.UpdateCursor намного быстрее, чем arcpy.CalculateField_management. Я написал скрипт, который работает на 55.000.000 конструктивных элементах, он был примерно в 5 раз медленнее с инструментом CalculateField.
offermann
Задача состоит в том, чтобы установить четыре подпроцесса и очистить память, так как это реальная точка сжатия. Как я обрисовал во втором абзаце, вы можете разделить подпроцессы по строкам, но это требует немного больше управления, чем один выбор.
MappaGnosis
1

Логика манипулирования данными может быть записана в виде оператора SQL UPDATE с использованием выражения CASE, которое можно выполнить с помощью GDAL / OGR, например, через OSGeo4W с gdal-filegdbустановленным.

Вот рабочий процесс, который использует osgeo.ogrвместо 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)

Для аналогичной таблицы, содержащей чуть более 1 миллиона записей, этот запрос занял 18 минут. Таким образом, обработка 16 миллионов записей может занять от 4 до 5 часов.

Майк Т
источник
К сожалению, скрипт является частью более крупного скрипта, написанного с использованием, arcpyно я ценю ответ. Я медленно пытаюсь использовать GDAL больше.
cptpython
1

Обновление кода в разделе № 2 по вашему вопросу не показывает, как вы присоединяете .csvфайл обратно к исходной таблице в вашей файловой базе геоданных. Вы говорите, что запуск вашего сценария занял ~ 5 минут. Это звучит справедливо, если вы только экспортировали .csvфайл без каких-либо соединений. Когда вы попытаетесь вернуть .csvфайл в ArcGIS, у вас возникнут проблемы с производительностью.

1) Вы не можете делать соединения напрямую с .csvтаблицей базы геоданных, потому что .csvфайл не имеет OID (поле, рассчитанное с уникальными значениями, не поможет, так как вам все равно потребуется преобразовать ваш .csvфайл в таблицу базы геоданных). Итак, несколько минут для Table To Tableинструмента GP (вы можете использоватьin_memory рабочее пространство для создания там временной таблицы, будет немного быстрее).

2) После загрузки .csvв таблицу базы геоданных вы захотите построить индекс для поля, в котором вы будете выполнять соединение (в вашем случае исходное значение objectidvaue из .csvфайла. Это займет несколько минут для таблицы с 16 млн строк).

3) Тогда вам нужно будет использовать Add Joinлибо Join Fieldинструменты GP, либо Ни один из них не будет хорошо работать на ваших больших столах.

4) После этого вам нужно использовать Calculate Fieldинструмент GP для вычисления вновь соединенных полей. Много минут здесь; Более того, вычисление полей занимает больше времени, когда поля, которые участвуют в вычислении, поступают из объединенной таблицы.

Одним словом, вы не получите ничего близкого к 5 минутам, которые вы упомянули. Если вы сделаете это через час, я был бы впечатлен.

Чтобы избежать обработки больших наборов данных в ArcGIS, я предлагаю перенести ваши данные вне ArcGIS во pandasфрейм данных и выполнить все ваши расчеты там. Когда вы закончите, просто запишите строки фрейма данных обратно в новую таблицу базы геоданных с помощью da.InsertCursor(или вы можете усечь существующую таблицу и записать свои строки в исходную).

Полный код, который я написал для тестирования, приведен ниже:

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

Ниже приведены выходные данные отладочного ввода-вывода (количество сообщений - это количество строк в используемой таблице) с информацией о времени выполнения для отдельных функций:

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

Вставка строки с помощью da.InsertCursorзанимает постоянное время, то есть если для вставки 1 строки требуется, скажем, 0,1 секунды, для вставки 100 строк потребуется 10 секунд. К сожалению, 95% + общего времени выполнения тратится на чтение таблицы базы геоданных и последующую вставку строк обратно в базу геоданных.

То же самое применимо к созданию pandasкадра данных из da.SearchCursorгенератора и к вычислению поля (полей). Поскольку число строк в вашей исходной таблице базы геоданных удваивается, увеличивается и время выполнения приведенного выше сценария. Конечно, вам все равно нужно использовать 64-битный Python, так как во время выполнения некоторые большие структуры данных будут обрабатываться в памяти.

Алекс Терешенков
источник
На самом деле, я собирался задать еще один вопрос, в котором говорилось бы об ограничениях метода, который я использовал, потому что я столкнулся с проблемами, которые вы рассмотрели выше, так что спасибо! Что я пытаюсь сделать: объединить четыре растра, а затем выполнить оператор if-else на основе столбцов, записать выходные данные в новый столбец и, наконец, выполнить Lookupсоздание растра на основе значений в новом столбце. В моем методе было много ненужных шагов и неэффективный рабочий процесс, я должен был упомянуть об этом в своем первоначальном вопросе. Живи и учись. Я попробую ваш сценарий позже на этой неделе.
cptpython