Usando o SRTM Global DEM para o cálculo da inclinação?

17

Eu baixei o SRTM GDEM (resolução de ~ 90 km).

Estou usando o ArcGIS 10.

Tentei usar o analista espacial para calcular a inclinação.

No entanto, não consigo calcular a inclinação.

Os valores de saída possuem apenas dois intervalos 0 e 0,1-90.

Não tenho muita certeza de qual é o problema?

user2543
fonte
Isso depende de onde você está analisando no mundo. Existem projeções diferentes para cada local. Onde você está examinando?
DJQ
5
A resolução é realmente ~ 90m, não ~ 90km.
Akheloes
Apenas um comentário, se você estiver em manutenção do Desktop, poderá fazer login no ArcGIS Online e usar seus serviços de elevação (sem necessidade de extensão de NA). A camada de inclinação é livre para usar como camada de referência. Na Austrália, temos os dados de 1 segundo do SRTM (~ 30m res) blogs.esri.com/esri/arcgis/2014/07/11/…
Simon

Respostas:

29

Esse parece ser um bom lugar para descrever uma maneira simples, rápida e mais do que razoavelmente precisa de calcular inclinações para um DEM globalmente extenso .

Princípios

Lembre-se de que a inclinação de uma superfície em um ponto é essencialmente a maior proporção de "elevação" para "rotação" encontrada em todos os rolamentos possíveis a partir desse ponto. O problema é que, quando uma projeção apresenta distorção de escala, os valores de "execução" serão computados incorretamente. Pior ainda, quando a distorção da balança varia com o rolamento - como é o caso de todas as projeções que não são conformes - como a inclinação varia com o rolamento será estimada incorretamente, impedindo a identificação precisa da relação máxima de subida: cálculo do aspecto).

Podemos resolver isso usando uma projeção conforme para garantir que a distorção da escala não varie com o rumo e, em seguida, corrigindo as estimativas de inclinação para levar em consideração a distorção da escala (que varia de ponto a ponto no mapa). O truque é usar uma projeção conforme global que permita uma expressão simples para sua distorção de escala.

A projeção de Mercator se ajusta à conta: assumindo que a escala está correta no Equador, sua distorção é igual à secante da latitude. Ou seja, as distâncias no mapa parecem ser multiplicadas pelo secante. Isso faz com que qualquer cálculo de inclinação calcule a subida: (seg (f) * corrida) (que é uma razão), onde f é a latitude. Para corrigir isso, precisamos multiplicar as inclinações computadas por segundo (f); ou, equivalentemente, divida-os por cos (f). Isso nos dá a receita simples:

Calcule a inclinação (como aumento: corrida ou uma porcentagem) usando uma projeção de Mercator, depois divida o resultado pelo cosseno da latitude.

Fluxo de trabalho

Para fazer isso com uma grade dada em graus decimais (como um SRTM DEM), execute as seguintes etapas:

  1. Crie uma grade de latitude. (Esta é apenas a grade de coordenadas y).

  2. Calcule seu cosseno.

  3. Projecto tanto o DM e o co-seno da latitude usando uma projecção Mercator na qual a escala é verdadeiro no equador.

  4. Se necessário, converta as unidades de elevação para concordar com as unidades das coordenadas projetadas (geralmente metros).

  5. Calcule a inclinação do DEM projetado como uma inclinação pura ou uma porcentagem ( não como um ângulo).

  6. Divida essa inclinação pela grade projetada de cosseno (latitude).

  7. Se desejar, reprojete a grade da inclinação em qualquer outro sistema de coordenadas para análise ou mapeamento adicional.

Os erros nos cálculos de inclinação serão de até 0,3% (porque esse procedimento usa um modelo de terra esférico em vez de um elipsóide, que é achatado em 0,3%). Esse erro é substancialmente menor que outros erros que entram nos cálculos de inclinação e, portanto, podem ser negligenciados.


Cálculos totalmente globais

A projeção Mercator não pode lidar com nenhum dos pólos. Para trabalhos em regiões polares, considere usar uma projeção estereográfica polar com escala real no polo. A distorção da escala é igual a 2 / (1 + sin (f)). Use esta expressão no lugar do item (f) no fluxo de trabalho. Especificamente, em vez de calcular uma grade de cosseno (latitude), calcule uma grade cujos valores sejam (1 + sin (latitude)) / 2 ( edit : use -latitude para o Polo Sul, conforme discutido nos comentários). Em seguida, proceda exatamente como antes.

Para uma solução global completa, considere dividir a grade terrestre em três partes - uma em torno de cada polo e outra em torno do equador -, executando um cálculo de inclinação separadamente em cada parte, usando uma projeção adequada, e fazendo um mosaico dos resultados. Um local razoável para dividir o globo é ao longo de círculos de latitude em latitudes de 2 * ArcTan (1/3), que é de cerca de 37 graus, porque nessas latitudes os fatores de correção Mercator e estereográficos são iguais entre si (tendo um valor comum de 5/4) e seria bom minimizar o tamanho das correções feitas. Como uma verificação dos cálculos, as grades devem estar em acordo muito próximo onde se sobrepõem (pequenas quantidades de imprecisão de ponto flutuante e diferenças devido à reamostragem das grades projetadas devem ser as únicas fontes de discrepâncias).

Referências

John P. Snyder, Projeções de Mapas - Um Manual de Trabalho . USGS Professional Paper 1395, 1987.

whuber
fonte
7
Encontro-me na posição, como tantas vezes faço, de agradecer mais uma vez a whuber por descrever uma solução e apresentar o raciocínio que a constrói. Eu tiro meu chapéu para você, senhor.
Matt Wilkie
Obrigado @matt. Não pretendi sugerir anteriormente que sua resposta (agora excluída) deve ser revogada: De fato, eu a votei porque você compartilhou um link para uma referência interessante do USGS que poderia ser útil para muitos leitores. (Meu comentário era crítica apenas uma passagem secundária em que o papel, não o próprio papel.)
whuber
ahh obrigado pelo esclarecimento. Eu já restaurou a resposta, confiando as pessoas têm informações suficientes na frente deles agora para fazer uma escolha informada :)
wilkie mate
2
Vindo de francês, demorei um pouco para traduzir a terminologia necessária para entender melhor essa ótima resposta, então achei que largar este link é uma boa ajuda para iniciantes como eu: webhelp.esri.com/arcgisdesktop/9.2/…
Akheloes
É uma ótima abordagem e eu já usei sua solução para gerar uma varredura de inclinação global. Uma dica da experiência prática: Como os valores de latitude sul do equador são negativos, você tem que usar o valor latitude absoluta na seguinte equação: (1 + sin (latitude)) / 2
Saleika
18

Resposta original

Suponho que as unidades horizontais da sua varredura estejam em graus ou em segundos. Você precisa reprojetar essa varredura em uma projeção espacial em que suas unidades horizontais e verticais são as mesmas (ou seja, se as unidades verticais estiverem em metros, sugiro usar o UTM, que possui unidades horizontais de metros).

Para reprojetar uma varredura com ArcCatalog / ArcGIS, consulte:

ArcToolbox> Ferramentas de Gerenciamento de Dados> Projeções e Transformações> Raster> Project Raster

Escolha uma referência espacial projetada que cubra sua região de interesse, por exemplo, tente uma zona UTM. Existem muitas outras opções que são melhor documentadas no manual . Observe que você não pode criar um conjunto de dados de declive para toda a Terra (se é isso que você está tentando fazer).

Melhor resposta, usando GDAL em uma escala

Agora que os dados do SRTM estão disponíveis globalmente , posso ver e trabalhar com os arquivos. O gdaldemutilitário da GDAL pode calcular declives e colinas usando uma opção de escala para uma proporção de unidades verticais para horizontais. O manual recomenda 111120 m / ° para algo como ladrilhos SRTM. Por exemplo, a partir de um shell OSGeo4W:

$ gdaldem slope -s 111120 -compute_edges N44E007.hgt N44E007_slope.tif

A -compute_edgesopção torna as bordas mais uniformes, se você deseja unir algumas peças. Ou calcule blocos para uma região grande. A desvantagem da técnica de "escala" é que as distâncias nas direções EW e NS não são iguais, exceto no equador, portanto, para os ladrilhos mais próximos dos pólos, pode haver algumas deturpações estranhas de inclinação.

Mike T
fonte
Vale ressaltar o seu comentário final: essa é uma solução ruim para pontos fora do Equador. Isso não é uma questão pequena de "deturpações estranhas": os resultados serão grosseiramente errados, especialmente em locais mais próximos dos poloneses do que do equador. A documentação para os gdaldemestados "Para locais fora do equador, seria melhor reprojetar sua grade usando gdalwarp antes de usar gdaldem." Infelizmente, isso não funcionará para conjuntos de dados que cobrem o globo, a menos que você os divida em pedaços pequenos (74 zonas UTM, talvez?), Projete-os, calcule as inclinações e faça mosaicos nos resultados.
whuber
7

Simplificando, não há um. Por definição, um sistema de coordenadas baseado em graus não é projetado. Na linguagem comum, dizemos que WGS84 é uma projeção "geográfica", mas isso é falso, apenas por conveniência.

Acho que me lembro de ler sobre um software ou processo para trabalhar com precisão com modelos de elevação no espaço geográfico não projetado, mas não consigo localizá-lo agora. De qualquer forma, seria experimental ou construa você mesmo a partir do tipo de processo do código.


Ahhh, encontrou: Desenvolvimento de um conjunto de dados de declive global para estimativa de ocorrência de deslizamentos de terra resultantes de terremotos (USGS). A página 4 descreve bem o problema

... o comprimento de um grau varia dependendo da sua localização latitudinal. No equador, um bloco de um grau por um grau é razoavelmente quadrado quando convertido em unidades de metros (111.321 metros na direção x por 110.567 metros na direção y ... mas mais próximo dos polos as distâncias na A direção x fica menor em função do cosseno da latitude, devido à convergência dos meridianos.A maioria dos pacotes GIS, incluindo o ArcGIS, opera apenas em pixels quadrados, usando um fator para ajustar as dimensões x, y ou z para uma unidade comum não é possível.

O artigo continua descrevendo os cálculos e ferramentas de software específicos ( , , ) que eles usaram para solucionar esse problema fundamental. O documento não inclui o código, mas, se solicitado com bom gosto, eles podem compartilhar. De qualquer forma, eu provavelmente perguntaria onde estão os resultados, sendo o USGS provavelmente já está online em algum lugar. :)

Matt Wilson
fonte
1
A sugestão desse artigo de que uma projeção equidistante azimutal poderia ser usada para calcular inclinações é equivocada e errada. De fato, dará declives corretos perto da origem da projeção, mas eles também ficarão progressivamente menos precisos à medida que a distância à origem aumenta.
whuber
obrigado por apontar isso. Os leitores, por favor, não deixe de ler gis.stackexchange.com/a/40464/108 , bem como, para o equilíbrio
Matt Wilkie
2

Parâmetros globais de DEM (onde a maioria das fórmulas são baseadas na suposição de espaço euclidiano) podem ser derivados com eficiência usando o sistema EQUI7 GRID (Bauer-Marschallinger et al. 2014). O EQUI7 GRID divide o mundo em 7 áreas terrestres, todas projetadas em um sistema de projeção equidistante com uma perda mínima de precisão. Veja um exemplo de DEM global com resolução de 250 m no EQUI7 GRID. Aqui você pode encontrar algum código de exemplo que mostra como derivar parâmetros globais de DEM usando SAGA GIS. Depois de concluir a derivação dos parâmetros DEM no sistema EQUI7 GRID, você pode transformar novamente todos os mapas em longlatcoordenadas WGS84 e criar um mosaico global usando GDAL.

Tom Hengl
fonte
Você poderia explicar como isso responde à pergunta? Se você está propondo usar projeções equidistantes para cálculos de inclinação, lembre-se de que é uma solução ruim devido às grandes distorções métricas relativas que ocorrem quando uma pessoa se afasta do centro da projeção. Embora focar sete dessas projeções em massas de terra ajude a aliviar esse problema, ainda não é a melhor escolha.
whuber
O artigo de Bauer-Marschallinger et al. (2014) explica por que essas projeções foram escolhidas para representar massas globais de terra (presume-se que elas tenham perda mínima de precisão). Concordo que qualquer projeção 2D acabará por levar a deformações, mas, tanto quanto eu sei, o EQUI7 é um bom compromisso entre perda de precisão e comodidade (álgebra 2D). Dito isto, hexágonos estão novamente sendo usados ​​para representar superfícies terrestres globais (embora a análise DEM com hexágonos 3D ainda seja um incômodo).
21136 Tom Hengl
Obrigado pela referência. Seu resumo sugere que ele resolve um problema bem diferente, o de "minimizar a amostragem excessiva de dados locais que aparecem durante a projeção de imagens de satélite genéricas em uma grade de varredura regular". Isso não implica que terá um bom desempenho para outros fins, como estimar declives.
whuber
É claro que o EQUI7 não resolve absolutamente o problema de estimar as inclinações locais com precisão, mas é provavelmente uma solução mais elegante do que usar a projeção de Mercator sugerida acima. Por fim, se desejar estimar declives com precisão perfeita, as únicas opções são provavelmente (1) usar projeções locais (equidistantes) por ladrilhos de menor tamanho (por exemplo, 100 por 100 km) com sobreposição de 10 a 20%, como mencionado também em Verdin et al. (2007) ou (2) para usar a grade hexagonal ( pacote dggridR ).
Tom Hengl
O problema não é precisão - está na produção de inclinações e aspectos sistematicamente tendenciosos. Como as projeções equidistantes distorcem diferencialmente as direções ortogonais às geodésicas originárias de seus centros, os aspectos sempre estarão errados (embora razoavelmente precisos perto dos centros onde toda a distorção é baixa) e os erros nas encostas crescerão rapidamente com a distância. É claro que o uso de muitas projeções locais funcionará, mas é exatamente o oposto da elegância que você valoriza.
whuber
-2

A inclinação é subir / correr. Calcule aumento e cálculo e você terá a sua resposta. É simples calcular a distância entre coordenadas geográficas. Isso apresentará menos erros de reamostragem em comparação à conversão para UTM, etc.

THK
fonte