Подсчет последовательных значений пикселей для набора растров с помощью ArcGIS Spatial Analyst?

23

Я использую ArcGIS 10 с Spatial Analyst.

У меня есть набор растров (всего 8), каждый из которых содержит 1 или 0 для каждой ячейки. Каждый растр представляет данные за разные годы. Для аргументов ради года 1 до года 8.

Я могу добавить все растры вместе, что даст мне окончательную сетку со значениями в диапазоне от 0 до 8. 8, указывающее, что ячейка была постоянно 1 для набора растров (все годы).

Я хотел бы найти для каждой ячейки самое длинное число последовательных единиц.

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

Мои навыки растровой обработки не так высоки, как мои навыки векторной обработки, и я хорошо посмотрел файл справки ESRI, но я не могу понять, как этого добиться, используя готовые инструменты геообработки?

Любые идеи?

Hornbydd
источник
1
Это на самом деле довольно крутой анализ. Как обычно, есть несколько способов сделать то, что вы пытаетесь сделать. Я думаю, что вам нужно будет немного программировать, чтобы пройтись по всем комбинациям.
MLowry
1
Общий комментарий (вдохновленный этим замечанием @MLowry): пожалуйста, задавайте вопросы, когда они выглядят интересными или четко сформулированы. Хорошие вопросы управляют всем на нашем сайте; пожалуйста, давайте сделаем все возможное, чтобы вознаградить тех, кто их спрашивает!
whuber

Ответы:

14

Поскольку это локальная операция, давайте разберемся, как это сделать для отдельной ячейки: Map Algebra позаботится обо всем остальном.

Во-первых, обратите внимание, что порядок растров имеет значение. Поэтому статистика единичных ячеек, такая как сумма ячеек, этого не сделает.

Если бы мы встретили последовательность, такую ​​как 01110101 в данной ячейке, мы бы обработали это от начала до конца и

  1. Начните с нуля.

  2. Увеличивайте счет каждый раз, когда мы сталкиваемся с 1.

  3. Сбрасывайте счетчик каждый раз, когда мы сталкиваемся с 0, после сохранения последнего счетчика .

  4. В конце возьмите максимальный сохраненный счет (включая окончательный счет).

Шаг 1 реализован с сеткой с постоянным нулем. Шаги 2 и 3 зависят от того, с чем мы сталкиваемся: поэтому это условная операция. Шаг 4 явно является локальным максимумом. Тогда мы бы закодировали это немного более формально:

count = 0
result = 0
For each value:
    If (value==1):
        count=count+1
    else
        result = max(result, count)
        count=0
result = max(result, count)

Лучше всего это сделать со скриптом Python, когда у вас много сеток, но с восьмью нет ничего сложного в том, чтобы развернуть цикл и выписать шаги вручную. Это обнаруживает небольшую проблему: result=max(longest,count)это своего рода «побочный эффект», который трудно кодировать растровыми операциями. (Но это можно сделать, как показано во втором решении ниже.) Это также неэффективно, поскольку добавляет дополнительные вычисления на каждом шаге. Поэтому мы немного изменим подход с целью отложить maxоперацию до конца. Это потребует сохранения отдельного счета на каждом этапе.

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

result1 = "grid1"
result2 = con("grid2"==1, "result1"+1, 0)
result3 = con("grid3"==1, "result2"+1, 0)
result4 = con("grid4"==1, "result3"+1, 0)
result5 = con("grid5"==1, "result4"+1, 0)
result6 = con("grid6"==1, "result5"+1, 0)
result7 = con("grid7"==1, "result6"+1, 0)
result8 = con("grid8"==1, "result7"+1, 0)
CellStatistics(["result1", "result2", "result3", "result4", "result5", "result6", "result7" "result8"], "max")

Фактический синтаксис зависит от вашей версии ArcMap. (Например, я CellStatisticsновичок в версии 10, но локальная максимальная операция всегда была доступна.)

В примере с вводом 01110101 последовательность сеток «result *» будет содержать значения 0, 1, 2, 3, 0, 1, 0, 1, поэтому в конце CellStatisticsбудет возвращаться значение 3, длина самой длинной строки 1 - х.

Если оперативной памяти недостаточно, решение может быть модифицировано для повторного использования промежуточных результатов, что приблизительно удваивает время выполнения:

result = "grid1"
temp = con("grid2"==1, "result"+1, 0)
result = CellStatistics[["temp", "result"], "max"]
temp = con("grid3"==1, "temp"+1, 0)
result = CellStatistics[["temp", "result"], "max"]
...
temp = con("grid8"==1, "temp"+1, 0)
CellStatistics[["temp", "result"], "max"]

В примере с вводом 01110101 значения ("temp", "result") будут (NoData, 0) после первой строки и после каждой пары операций ("Con", "CellStatistics") значения будут равны (1 , 1), (2, 2), (3, 3), (0, 3), (1, 3), (0, 3), (1, 3). Еще раз окончательное значение 3.

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

Whuber
источник
В блоке кода типа может быть опечатка: count = count = 1, вероятно, должно быть count = count + 1
MLowry
1
@ML Спасибо (хорошие глаза!): Это исправлено. Трудно сделать псевдокод абсолютно правильным; человеческий обзор является реальным активом в поиске ошибок. Кроме того, хотя я не тестировал решения в ArcGIS, я реализовал первое решение в R, поэтому у меня есть уверенность, что подход правильный.
whuber
1
"еще раз", ты снова тот человек, который знает! Боже, помоги остальным из нас, если тебя когда-нибудь переедет автобус! Ваш первоначальный подход к Python был тем направлением, о котором я думал, но я знаю, что с растрами часто можно сделать все настолько гладко, как вы это доказали. Если вы окажетесь в Британии, будет честью купить вам пинту нашего лучшего плоского пива комнатной температуры! :)
Хорнбидд
Спасибо, Дункан: но зацени отличное решение Энди Харфута!
whuber
14

Просто поговорите об этом и спросите себя, можете ли вы подойти к проблеме, рассматривая входные сетки как двоичный поток. Это позволит вам объединить их, чтобы получить уникальное итоговое целое число для последовательности - т.е. 01110101 = 117. Затем это значение можно было бы реклассифицировать, чтобы получить максимальное количество последовательных единиц.

Вот пример, показывающий один из способов объединения восьми сеток:

2*(2*(2*(2*(2*(2*(2*"g8" + "g7") + "g6") + "g5") + "g4") + "g3") + "g2") + "g1"

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

В таблице переклассификации должны быть указаны максимальные длины прогонов для всех значений от 00000000B = 0 до 11111111B = 255. По порядку они следующие:

0, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 6, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 4, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 7, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 6, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 7, 8

Этот подход ограничен примерно 20 сетками в ArcGIS: использование большего количества может создать громоздкую таблицу атрибутов. ( Combineспециально ограничено 20 сетками.)

Энди Харфут
источник
Огромный +1: это очень хорошая идея. (Единственное ограничение заключается в том, что когда задействовано более 31 сетки, у вас не хватит битов для использования.) Я позволил себе немного конкретизировать вашу идею, чтобы другие могли увидеть, насколько легко ее реализовать.
whuber
3

Задумывались ли вы об изменении значений от 0 и 1 до значений со степенью 2 (1,2,4,8,16,32). Когда вы объедините 8 сеток, вы получите уникальные значения для каждой ячейки, которые предоставят вам последовательную информацию (то есть: значение 3 означает год 1 и 2, где значение 54 будет годами 6-8).

Просто мысль

Райан Гарнетт
источник
Это именно то, что @ Энди Харфут предложил несколько часов назад, Райан. :-)
whuber
Спасибо и извините. Я прочитал это на моем телефоне в отпуске.
Райан Гарнетт