Используя ArcGIS 10, у меня есть растр, в котором я хотел бы найти пиксель с максимальным значением в растре и вернуть его местоположение (центр пикселя) в десятичных градусах. Я хотел бы повторить этот процесс, возвращая местоположение второго по величине значения растра, затем третьего и т. Д., Чтобы в конце я получил список из N местоположений, которые имеют самые высокие значения в растре по порядку.
Я полагаю, что это легче всего сделать с помощью скрипта Python, но я открыт для других идей, если есть лучший способ.
Ответы:
Если вы счастливы использовать R , есть пакет под названием растр . Вы можете прочитать в растре, используя следующую команду:
Затем, когда вы идете, чтобы посмотреть на него (набрав
test
), вы можете увидеть следующую информацию:Могут быть и более эффективные способы манипулирования растром, но одним из способов поиска нужной информации может быть нахождение наибольшего значения и определение его местоположения в матрице, а затем добавление этого к нижним экстентам.
источник
R
, вы можете использовать стандартныеR
функции илиgetValues
метод для доступа к значениям ячеек. Отсюда легко определить самые высокие значения и их местоположение.Ответ может быть получен путем объединения индикаторной сетки верхних 1% значений с сетками широты и долготы. Хитрость заключается в создании этой сетки индикаторов, потому что ArcGIS (все же! После 40 лет!) Не имеет процедуры для ранжирования растровых данных.
Одно из решений для растров с плавающей точкой является итеративным, но, к счастью, быстрым . Пусть n будет количеством ячеек данных. Эмпирическое распределение значений состоит из всех пар (г, п (г)) , где г представляет собой значение в сетке и п (г) представляет собой количество ячеек в сетке со значениями меньше или равен г . Мы получаем кривую, соединяющую (-infinity, 0) с (+ infinity, n) из последовательности этих вершин, упорядоченной по z . Тем самым он определяет функцию f , где (z, f (z)) всегда лежит на кривой. Вы хотите найти точку (z0, 0,99 * n) на этой кривой.
Другими словами, задача состоит в том, чтобы найти ноль f (z) - (1-0.01) * n . Сделайте это с любой процедурой поиска нуля (которая может обрабатывать произвольные функции: эта не дифференцируема). Самым простым, который часто эффективен, является угадывание и проверка: изначально вы знаете, что z0 лежит между минимальным значением zMin и максимальным zMax. Угадайте любую разумную ценность строго между этими двумя. Если предположение слишком низкое, установите zMin = z0; в противном случае установите zMax = z0. Теперь повторите. Вы быстро дойдете до решения; вы достаточно близко, когда zMax и zMin достаточно близко. Чтобы быть консервативным, выберите окончательное значение zMin в качестве решения: оно может набрать несколько дополнительных очков, которые вы можете отбросить позже. Для более сложных подходов, см. Главу 9 Числовых Рецептов (ссылка идет на старую бесплатную версию).
Оглядываясь назад на этот алгоритм, вы обнаруживаете, что вам необходимо выполнять только два вида растровых операций : (1) выбрать все ячейки, меньшие или равные некоторому целевому значению, и (2) подсчитать выбранные ячейки. Это одни из самых простых и быстрых операций. (2) может быть получено в виде зонального счета или путем считывания одной записи из таблицы атрибутов сетки выбора.
источник
Я сделал это некоторое время назад, хотя мое решение использует GDAL (так что это не только для ArcGIS). Я думаю, что вы можете получить массив NumPy из растра в ArcGIS 10, но я точно не знаю. NumPy обеспечивает простую и мощную индексацию массивов, как
argsort
и другие. Этот пример не обрабатывает NODATA и не преобразует координаты из проецируемой в широту / долготу (но это не сложно сделать с помощью osgeo.osr, предоставляемого с GDAL)Показывает следующее для моего тестового растрового файла:
источник
NODATA = rast_band.GetNoDataValue()
, затем с использованием либо значения NaN (rast[rast == NODATA] = np.nan
), либо с помощью маскированного массива (rast = np.ma.array(rast, mask=(rast == NODATA))
). Более сложный трюк состоит в том,argsort
чтобы каким-то образом удалить значения NODATA из анализа или просто пропустить их в цикле for, если они являются NaN / замаскированными.