Поиск местоположений с самыми высокими значениями в растре с помощью ArcGIS Desktop?

12

Используя ArcGIS 10, у меня есть растр, в котором я хотел бы найти пиксель с максимальным значением в растре и вернуть его местоположение (центр пикселя) в десятичных градусах. Я хотел бы повторить этот процесс, возвращая местоположение второго по величине значения растра, затем третьего и т. Д., Чтобы в конце я получил список из N местоположений, которые имеют самые высокие значения в растре по порядку.

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

мГА
источник
Вы пытались преобразовать сетку в точки, а затем добавить поля X, Y и сортировку?
Якуб Сисак GeoGraphics
Растровые значения являются числами с плавающей точкой или целыми числами?
whuber
@ Якуб - Нет, не знаю. Возможно, меня заинтересуют только первые 1% или около того точек, поэтому я не знаю, стоит ли добавлять поля x, y для всех точек, а затем сортировать. Может быть, если нет более эффективного варианта?
MGA
@whuber - растровые значения являются числами с плавающей точкой.
MGA,
@ Мга, стоит попробовать. Преобразование происходит довольно быстро, и добавление XY также является инструментом по умолчанию. Удаление ненужных записей - прямая операция, и все они могут быть объединены в одну модель. Просто идея.
Якуб Сисак GeoGraphics

Ответы:

5

Если вы счастливы использовать R , есть пакет под названием растр . Вы можете прочитать в растре, используя следующую команду:

install.packages('raster')
library(raster)
test <- raster('F:/myraster')

Затем, когда вы идете, чтобы посмотреть на него (набрав test), вы можете увидеть следующую информацию:

class       : RasterLayer 
dimensions  : 494, 427, 210938  (nrow, ncol, ncell)
resolution  : 200, 200  (x, y)
extent      : 1022155, 1107555, 1220237, 1319037  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
values      : F:/myraster 
min value   : 0 
max value   : 1 

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

djq
источник
1
+1 За эту ссылку. После того как вы прочитали растр R, вы можете использовать стандартные Rфункции или getValuesметод для доступа к значениям ячеек. Отсюда легко определить самые высокие значения и их местоположение.
whuber
1
Благодаря вашей рекомендации, это то, что я в итоге сделал. Использование растрового пакета в R было проще простого по сравнению с этим в ArcGIS. Я также продолжил использовать другие пространственные анализы в R и был очень доволен результатами. Отличный совет!
MGA
8

Ответ может быть получен путем объединения индикаторной сетки верхних 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) может быть получено в виде зонального счета или путем считывания одной записи из таблицы атрибутов сетки выбора.

Whuber
источник
7

Я сделал это некоторое время назад, хотя мое решение использует GDAL (так что это не только для ArcGIS). Я думаю, что вы можете получить массив NumPy из растра в ArcGIS 10, но я точно не знаю. NumPy обеспечивает простую и мощную индексацию массивов, как argsortи другие. Этот пример не обрабатывает NODATA и не преобразует координаты из проецируемой в широту / долготу (но это не сложно сделать с помощью osgeo.osr, предоставляемого с GDAL)

import numpy as np
from osgeo import gdal

# Open raster file, and get GeoTransform
rast_src = gdal.Open(rast_fname)
rast_gt = rast_src.GetGeoTransform()

def get_xy(r, c):
    '''Get (x, y) raster centre coordinate at row, column'''
    x0, dx, rx, y0, ry, dy = rast_gt
    return(x0 + r*dx + dx/2.0, y0 + c*dy + dy/2.0)

# Get first raster band
rast_band = rast_src.GetRasterBand(1)

# Retrieve as NumPy array to do the serious work
rast = rast_band.ReadAsArray()

# Sort raster pixels from highest to lowest
sorted_ind = rast.argsort(axis=None)[::-1]

# Show highest top 10 values
for ind in sorted_ind[:10]:
    # Get row, column for index
    r, c = np.unravel_index(ind, rast.shape)
    # Get [projected] X and Y coordinates
    x, y = get_xy(r, c)
    print('[%3i, %3i] (%.3f, %.3f) = %.3f'%
          (r, c, x, y, rast[r, c]))

Показывает следующее для моего тестового растрового файла:

[467, 169] (2813700.000, 6353100.000) = 844.538
[467, 168] (2813700.000, 6353200.000) = 841.067
[469, 168] (2813900.000, 6353200.000) = 840.705
[468, 168] (2813800.000, 6353200.000) = 840.192
[470, 167] (2814000.000, 6353300.000) = 837.063
[468, 169] (2813800.000, 6353100.000) = 837.063
[482, 166] (2815200.000, 6353400.000) = 833.038
[469, 167] (2813900.000, 6353300.000) = 832.825
[451, 181] (2812100.000, 6351900.000) = 828.064
[469, 169] (2813900.000, 6353100.000) = 827.514
Майк Т
источник
+1 Спасибо, что поделились этим. Я не вижу невозможности обрабатывать NoData как ограничение: просто преобразуйте все NoData в крайне отрицательные значения, прежде чем продолжить. Также обратите внимание, что если произойдет какое-либо перепроецирование сетки , ответы, скорее всего, изменятся из-за повторной выборки сетки, поэтому обычно не требуется, чтобы перепроецирование происходило автоматически во время подобных вычислений. Вместо этого сообщенные координаты могут быть впоследствии перепроектированы. Таким образом, ваше решение является совершенно общим.
whuber
Обработка NODATA может быть реализована сначала путем получения значения из растра NODATA = rast_band.GetNoDataValue(), затем с использованием либо значения NaN ( rast[rast == NODATA] = np.nan), либо с помощью маскированного массива ( rast = np.ma.array(rast, mask=(rast == NODATA))). Более сложный трюк состоит в том, argsortчтобы каким-то образом удалить значения NODATA из анализа или просто пропустить их в цикле for, если они являются NaN / замаскированными.
Майк Т