Encontrando locais de valores mais altos em varredura usando o ArcGIS Desktop?

12

Usando o ArcGIS 10, tenho uma varredura na qual gostaria de encontrar o pixel com o valor máximo na varredura e retornar sua localização (centro do pixel) em graus decimais. Gostaria de percorrer esse processo retornando a localização do segundo valor mais alto da varredura, depois do terceiro e assim por diante, para que, no final, eu tenha uma lista de N localizações que tenham os valores mais altos na varredura em ordem.

Imagino que isso possa ser feito com mais facilidade usando um script Python, mas estou aberto a outras idéias, se houver uma maneira melhor.

mga
fonte
você tentou converter a grade em pontos e adicionar campos X, Y e classificação?
Jakub Sisak GeoGraphics
Os valores da varredura são flutuantes ou inteiros?
whuber
@ Jakub - Não, eu não tenho. Provavelmente, só vou me interessar entre os 1% ou mais dos pontos, por isso não sei se vale a pena adicionar campos x, y para todos os pontos e depois classificá-los. Talvez se não houver uma opção mais eficiente?
mga
@whuber - Os valores da varredura são flutuadores.
mga
@ga, vale a pena tentar. A conversão é bastante rápida e adicionar o XY também é uma ferramenta padrão. A exclusão de registros indesejados é uma operação direta e todos podem ser agrupados em um único modelo. Apenas uma ideia.
Jakub Sisak GeoGraphics

Respostas:

5

Se você gosta de usar o R , existe um pacote chamado raster . Você pode ler em uma varredura usando o seguinte comando:

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

Então, quando você olhar (digitando test), poderá ver as seguintes informações:

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 

Pode haver maneiras melhores de manipular a varredura, mas uma maneira de encontrar as informações desejadas pode ser encontrar o valor mais alto e obter a localização da matriz, e depois adicioná-lo às extensões mais baixas.

djq
fonte
1
+1 para esta referência. Depois de ler a varredura R, você pode usar Rfunções padrão ou o getValuesmétodo para acessar os valores das células. A partir daí, é fácil identificar os valores mais altos e suas localizações.
whuber
1
Graças à sua recomendação, é isso que acabei fazendo. Usar o pacote raster no R foi muito fácil comparado a tentar isso no ArcGIS. Também continuei usando outras análises espaciais em R e fiquei muito feliz com os resultados. Ótimo conselho!
mga
8

A resposta pode ser obtido por combinação de uma grelha de indicador de topo de 1% de valores com grades de latitude e longitude. O truque está na criação dessa grade de indicadores, porque o ArcGIS (ainda! Após 40 anos!) Não possui um procedimento para classificar dados rasterizados.

Uma solução para rasters de ponto flutuante é iterativa, mas felizmente rápida . Seja n o número de células de dados. A distribuição cumulativa empírica dos valores consiste em todos os pares (z, n (z)) em que z é um valor na grade en (z) é o número de células na grade com valores menores ou iguais a z . Temos uma curva conectando (-in infinito, 0) a (+ infinito, n) a partir da sequência desses vértices ordenada por z . Assim, define uma função f , onde (z, f (z)) sempre fica na curva. Você deseja encontrar um ponto (z0, 0,99 * n) nesta curva.

Em outras palavras, a tarefa é encontrar um zero de f (z) - (1-0.01) * n . Faça isso com qualquer rotina de busca zero (que pode lidar com funções arbitrárias: essa não é diferenciável). O mais simples, que geralmente é eficiente, é adivinhar e verificar: inicialmente, você sabe que z0 fica entre o valor mínimo zMin e o máximo zMax. Adivinhe qualquer valor razoável estritamente entre esses dois. Se a estimativa for muito baixa, defina zMin = z0; caso contrário, defina zMax = z0. Agora repita. Você convergirá rapidamente para a solução; você está próximo o suficiente quando o zMax e o zMin estão suficientemente próximos. Para ser conservador, escolha o valor final de zMin como a solução: ele pode obter alguns pontos extras que você pode descartar mais tarde. Para abordagens mais sofisticadas, consulte o Capítulo 9 de Receitas Numéricas (o link vai para uma versão gratuita mais antiga).

Analisando esse algoritmo, você precisa executar apenas dois tipos de operações de varredura : (1) selecione todas as células menores ou iguais a algum valor-alvo e (2) conte as células selecionadas. Essas estão entre as operações mais simples e rápidas. (2) pode ser obtido como uma contagem zonal ou lendo um registro da tabela de atributos da grade de seleção.

whuber
fonte
7

Eu fiz isso há algum tempo, embora minha solução esteja usando GDAL (então, isso não é apenas para ArcGIS). Eu acho que você pode obter uma matriz NumPy de uma varredura no ArcGIS 10, mas não tenho certeza. O NumPy fornece indexação de matriz simples e poderosa, como argsortoutras. Este exemplo não lida com NODATA ou transforma coordenadas de projetadas em latas / longas (mas isso não é difícil de fazer com osgeo.osr, fornecido com 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]))

Mostra o seguinte para o meu arquivo raster de teste:

[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
Mike T
fonte
+1 Obrigado por compartilhar isso. Não vejo a incapacidade de tratar o NoData como uma limitação: basta converter todos os NoData em valores extremamente negativos antes de prosseguir. Observe também que, se ocorrer alguma reprojeção da grade , as respostas provavelmente mudarão devido à reamostragem da grade; portanto, normalmente não se deseja que a reprojeção ocorra automaticamente durante um cálculo como esse. Em vez disso, as coordenadas relatadas podem ser reprojetadas posteriormente. Portanto, sua solução é perfeitamente geral.
whuber
A manipulação do NODATA pode ser implementada obtendo o valor da varredura primeiro e NODATA = rast_band.GetNoDataValue(), em seguida, usando um valor NaN ( rast[rast == NODATA] = np.nan) ou usando uma matriz mascarada ( rast = np.ma.array(rast, mask=(rast == NODATA))). O truque mais complicado é conseguir, argsortde alguma forma, remover os valores NODATA da análise ou simplesmente ignorá-los no loop for se estiverem NaN / mascarados.
Mike T