Как реклассифицировать очень большой набор данных о земном покрове?

10

Рассмотрим набор данных NLCD2001 Land Cover для Аляски ( ссылка для скачивания ). Мне нужно переклассифицировать этот набор данных, чтобы сохранить только пиксели со значениями 41, 42 и 43; все остальные значения пикселей должны стать NoData (или 0, если необходимо).

Это кажется простой задачей, требующей всего одного вызова инструмента переклассификации. К сожалению, каждый вызов приводит к неопределенному и бесполезному сообщению об ошибке:

Executing: Reclassify "D:\ak_nlcd_2001_land_cover_3-13-08_se5.img" Value "0 40 0;41 41;42 42;43 43;44 255 0;NODATA 0" "D:\alaska_reclassified.tif" DATA 
Start Time: Thu Jan 03 09:23:13 2013
ERROR 999998: Unexpected Error.
Failed to execute (Reclassify).
Failed at Thu Jan 03 09:23:13 2013 (Elapsed Time: 0.00 seconds)

Как я могу реклассифицировать этот набор растровых данных? Я использую ArcCatalog 10.0, Build 4000, с включенным расширением Spatial Analyst.

DoggoDougal
источник
Извлечение по атрибутам, кажется, также делает то, что мне нужно, но, к сожалению, приводит к еще одной «непредвиденной ошибке».
DoggoDougal
Возможно, пробовал другой набор данных? Два процесса, терпящие неудачу в одном и том же наборе данных, заставляют задуматься ...
Чед Купер
2
Обычно это reclassifyдолжно быть последнее средство, потому что оно настолько общее по объему, что, вероятно, использует методы, которые менее эффективны, чем те, которые можно получить, когда реклассификацию легко выразить арифметически или логически. В данном случае критерий переклассификации настолько прост, что вам следует сначала попробовать его Conили даже выполнить арифметические операции (потому что они быстрые). Например, "grid" * ("grid" >= 41) * ("grid" <= 43)должен это сделать. ОЗУ не должно быть проблемой - Spatial Analyst автоматически запускает свой растровый ввод-вывод, и это локальные операции.
whuber
1
Inlistэто хорошее решение (+1). Я мог использовать conи контролировать использование оперативной памяти во время операции. Он никогда не превышал 180 МБ, что чуть больше, чем оперативная память, используемая только для запуска ArcMap. Плитка в ArcGIS автоматическая - вы даже не можете управлять ей (если вы не программируете интерфейс C / Fortran). Похоже, что ограничения ОЗУ не представляют большой проблемы.
whuber
1
@whuber, conу меня тоже работал, с условием "Value" >= 41 AND "Value" <= 43. Я бы пошел с этим решением, но я не уверен, будут ли дополнительные растровые значения представлять интерес в будущем. Очевидно, я мог бы добавить ORк предложению where, но тогда это стало бы более сложным. InListкажется наиболее простым решением в отношении удобочитаемости и удобства обслуживания.
DoggoDougal

Ответы:

9

Первый прикрепленный скрипт успешно переклассифицировал ваши данные AK NLCD примерно за 15 минут (i7, 12 ГБ ОЗУ). Поскольку исходный набор данных составляет почти 7 ГБ, вы можете столкнуться с проблемами с памятью. Если вы не можете обработать весь набор данных в одном фрагменте, попробуйте разделить его со вторым сценарием до переклассификации. Я рекомендую взять небольшое подмножество данных (щелкните растровый слой правой кнопкой мыши в TOC> Данные> Экспорт данных> Экстент (фрейм данных) и протестируйте первый сценарий. Как только вы наберете параметры для команды переклассификации, затем перейдите к переклассификации весь набор данных или его разделение. Или попробуйте загрузить 64-битный продукт фоновой геообработки для ArcGIS 10.1 SP1, доступный здесь . Удачи.

Сценарий 1

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Overwrite output
env.overwriteOutput = 1

# Set environment settings
env.workspace = r'C:\temp'
Dir = env.workspace

# Set local variables
inRaster = Dir + "\\" + "nlcd_subset.img"
reclassField = "Value"
remap = RemapValue([[0, 40, 0], [41, 41],[42,42], [43,43], [44, 256, 0]])

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute Reclassify
outReclassify = Reclassify(inRaster, reclassField, remap, "NODATA")

# Save the output 
outReclassify.save(r"C:\temp\nlcd_test.img")

Изменить : Если вам нужно разделить ваши данные перед обработкой, этот скрипт должен помочь:

Сценарий 2

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Overwrite output
env.overwriteOutput = 1

# Set environment settings
env.workspace = r'C:\temp'
Dir = env.workspace

# Set local variables
inRaster = Dir + "\\" "nlcd" + "\\" + "nlcd_ak.img"
outFolder = Dir
reclassField = "Value"
remap = RemapValue([[0, 40, 0], [41, 41],[42,42], [43,43], [44, 256, 0]])

# Split Rasters
# Equally split a large TIFF image by number of images
arcpy.SplitRaster_management(inRaster, outFolder, "split", "NUMBER_OF_TILES", "#",
                             "NEAREST", "2 2", "#", "4", "PIXELS",\
                             "#", "#")

# List rasters for processing
rasters = arcpy.ListRasters()


for ras in rasters:
    print "processing..." + ras

    # Define new name
    name = "class_" + ras  

    # Execute Reclassify
    outReclassify = Reclassify(ras, reclassField, remap, "NODATA")

    # Save the output 
    outReclassify.save(Dir + "\\" + name)
Аарон
источник
3
С точки зрения производительности, было бы интересно попробовать альтернативный подход с использованием arcpy.RasterToNumPyArray () и выполнить переклассификацию в numpy. Возможно, вы все равно захотите разбить растр на тайлы в целях памяти, но я знаю, что с GDAL переклассификация массивов-пустышек выполняется очень быстро.
DavidF
@DavidF Согласен, вероятно, будет значительное улучшение производительности.
Аарон
Спасибо за советы, Аарон. Я выполню его, как только закончу другой обходной путь, который, кажется, требует удаления карты цветов ( ссылка здесь ). Этот метод также требует разбиения растра, поэтому меня интересует, не произошла ли переклассификация оригинала из-за использования памяти или по какой-либо другой причине.
DoggoDougal
@torik Нет проблем - я с удовольствием отдаю свои два цента. Я думаю, что удаление карты цветов это не тот путь. Скорее, я бы сосредоточился на разделении данных или 64-битной фоновой обработке.
Аарон
@ Аарон, учитывая, что вы предоставили код для выполнения мозаики, как вы создали растр подмножества, который вы использовали для получения изображенных результатов? Я завершил листинг SplitRaster (создав 100 подмножеств всего набора растровых данных) и попытался перебрать их все для переклассификации. К сожалению, переклассификация не удалась, что привело к тому же сообщению «Неожиданная ошибка».
DoggoDougal
4

whuber прокомментировал использование логических инструментов для выражения этой реклассификации . Немного покопавшись , я обнаружил, что InList , как часть набора инструментов логической математики Spatial Analyst, удовлетворяет мои потребности.

import arcpy

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
from arcpy.sa import InList

# Pixel values of interest, named according to Table 2 of
#  http://landcover.usgs.gov/pdf/anderson.pdf
DECIDUOUS_FOREST = 41
EVERGREEN_FOREST = 42
MIXED_FOREST = 43

inRaster = r'D:\AK_NLCD_2001_land_cover_3-13-08\ak_nlcd_2001_land_cover_3-13-08_se5.img'
accepted_raster_values = [DECIDUOUS_FOREST, EVERGREEN_FOREST, MIXED_FOREST]
filteredAlaska = InList(inRaster, accepted_raster_values)
filteredAlaska.save(r'C:\alaska\ak_woods')

На сегодняшний день это самое простое решение, которое я смог найти, оно работает быстрее всего и не требует учета листов исходного набора данных. Нет необходимости учитывать доступную оперативную память машины, поскольку этот инструмент будет считывать данные прямо с диска и сохранять результаты прямо на диске.

Отфильтрованный результат Аляски с использованием InList

DoggoDougal
источник
+1 Молодец и отличное решение. Из любопытства, сколько времени заняла обработка?
Аарон
@ Аарон, обработка всей Аляски занимает 13 минут и 23,4 секунды. Образец подмножество , который является одним из 100 одинакового размера подмножеств , созданных SplitRaster_management, занимает 7,04 секунды.
DoggoDougal
Интересно, что примерно одинаковое время обработки между двумя методами (т. Е. Предполагается, что мы работали с одинаковыми системами).
Аарон
У меня Intel Core 2 Duo E6850 с тактовой частотой 3 ГГц, 4 ГБ оперативной памяти, работающий под управлением 64-разрядной версии Windows 7. Вскоре я сделаю анализ времени вашего решения. Я пока застрял в Arc 10.0, иначе я бы изучил 64-битную фоновую обработку.
DoggoDougal
1

Я использовал набор данных, упомянутый в исходном сообщении, с версией arcmap 10.4 dev. Переклассификация завершается неудачно, когда выходной растр является сеткой, поскольку количество переклассифицированных ячеек выходит за пределы того, что может быть сохранено в поле COUNT НДС сетки. Когда выходной растр представляет собой fgdb, он успешно выполняется для меня примерно за 11 минут на старом 4-ядерном компьютере под управлением Windows 8. Неформатные растровые форматы должны работать, так как они используют плавающие значения двойной точности для поля счетчика. Я ожидаю, что вы должны получить то же поведение с 10.2 или 10.3 выпущенными версиями. Мы рассмотрим использование другого растрового формата для вывода по умолчанию для Reclassify.

jt64
источник