Como reclassificar um conjunto de dados de cobertura do solo muito grande?

10

Considere o conjunto de dados NLCD2001 Land Cover para o Alasca ( link para download ). Preciso reclassificar esse conjunto de dados para que apenas os pixels dos valores 41, 42 e 43 sejam preservados; todos os outros valores de pixel devem se tornar NoData (ou 0, se necessário).

Parece uma tarefa simples, exigindo apenas uma chamada para a ferramenta Reclassificar. Infelizmente, todas as chamadas resultam em uma mensagem de erro vaga e inútil:

Executing: Reclassify "D:\ak_nlcd_2001_land_cover_3-13-08_se5.img" Value "0 40 0;41 41;42 42;43 43;44 255 0;NODATA 0" "D:\alaska_reclassified.tif" DATA 
Start Time: Thu Jan 03 09:23:13 2013
ERROR 999998: Unexpected Error.
Failed to execute (Reclassify).
Failed at Thu Jan 03 09:23:13 2013 (Elapsed Time: 0.00 seconds)

Como posso reclassificar esse conjunto de dados raster? Estou usando o ArcCatalog 10.0, Build 4000, com a extensão Spatial Analyst ativada.

DoggoDougal
fonte
Extrair por atributos parece também fazer o que eu preciso, mas infelizmente resulta em outro "Erro inesperado".
precisa saber é o seguinte
Tentei outro conjunto de dados, talvez? Dois processos falhando no mesmo conjunto de dados torna ya saber ...
Chad Cooper
2
Normalmente, reclassifydeve ser um último recurso, porque é de escopo tão geral que provavelmente usa métodos menos eficientes do que os que podem ser obtidos quando a reclassificação é fácil de expressar aritmeticamente ou logicamente. No presente caso, o critério para a reclassificação é tão simples que você deve experimentá-lo primeiro Conou mesmo com operações aritméticas diretas (porque são rápidas). Por exemplo, "grid" * ("grid" >= 41) * ("grid" <= 43)deveria fazê-lo. A RAM não deve ser um problema - o Spatial Analyst abre automaticamente sua E / S raster e essas são operações locais.
whuber
1
Inlisté uma boa solução (+1). Consegui usar cone monitorar o uso da RAM durante a operação. Nunca excedeu 180 MB, pouco mais que a RAM usada apenas para iniciar o ArcMap. O mosaico no ArcGIS é automático - você nem consegue controlá-lo (a menos que esteja programando na interface C / Fortran). Parece que as limitações de RAM são de pouca preocupação.
whuber
1
@ Whuber, contrabalhou para mim também, com a condição "Value" >= 41 AND "Value" <= 43. Eu teria optado por essa solução, mas não tenho certeza se valores adicionais de varredura serão de interesse no futuro. Obviamente, eu poderia adicionar uma ORcláusula where, mas ela se tornaria mais complicada. InListparece a solução mais direta em relação à legibilidade e manutenção.
precisa saber é o seguinte

Respostas:

9

O primeiro script anexado reclassificou com êxito seus dados do AK NLCD em cerca de 15 minutos (i7, máquina de 12 GB de RAM). Como o conjunto de dados original tem quase 7 GB, você pode encontrar problemas de memória. Se você não puder processar o conjunto de dados inteiro em um pedaço, tente dividi-lo com o segundo script antes da reclassificação. Minha recomendação é pegar um pequeno subconjunto de dados (clique com o botão direito do mouse na camada raster em TOC> Dados> Exportar dados> Extensão (quadro de dados) e testar o primeiro script. Depois de discar os parâmetros para o comando reclassify, avance para reclassificar Como alternativa, tente fazer o download do produto de geoprocessamento em segundo plano de 64 bits para ArcGIS 10.1 SP1, disponível aqui . Boa sorte.

Script 1

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Overwrite output
env.overwriteOutput = 1

# Set environment settings
env.workspace = r'C:\temp'
Dir = env.workspace

# Set local variables
inRaster = Dir + "\\" + "nlcd_subset.img"
reclassField = "Value"
remap = RemapValue([[0, 40, 0], [41, 41],[42,42], [43,43], [44, 256, 0]])

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute Reclassify
outReclassify = Reclassify(inRaster, reclassField, remap, "NODATA")

# Save the output 
outReclassify.save(r"C:\temp\nlcd_test.img")

Editar : se você precisar dividir seus dados antes do processamento, este script deve ajudar:

Script 2

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Overwrite output
env.overwriteOutput = 1

# Set environment settings
env.workspace = r'C:\temp'
Dir = env.workspace

# Set local variables
inRaster = Dir + "\\" "nlcd" + "\\" + "nlcd_ak.img"
outFolder = Dir
reclassField = "Value"
remap = RemapValue([[0, 40, 0], [41, 41],[42,42], [43,43], [44, 256, 0]])

# Split Rasters
# Equally split a large TIFF image by number of images
arcpy.SplitRaster_management(inRaster, outFolder, "split", "NUMBER_OF_TILES", "#",
                             "NEAREST", "2 2", "#", "4", "PIXELS",\
                             "#", "#")

# List rasters for processing
rasters = arcpy.ListRasters()


for ras in rasters:
    print "processing..." + ras

    # Define new name
    name = "class_" + ras  

    # Execute Reclassify
    outReclassify = Reclassify(ras, reclassField, remap, "NODATA")

    # Save the output 
    outReclassify.save(Dir + "\\" + name)
Aaron
fonte
3
Do ponto de vista do desempenho, seria interessante tentar uma abordagem alternativa usando arcpy.RasterToNumPyArray () e fazer a reclasse em numpy. Você provavelmente desejaria dividir a varredura em blocos de qualquer maneira para fins de memória, mas eu sei que com o GDAL, a reclassificação de matrizes numpy é muito rápida.
DavidF 03/01
@DavidF Concordado, provavelmente haveria uma melhoria significativa no desempenho.
Aaron
Obrigado pelas dicas, Aaron. Vou executá-lo assim que terminar outra solução alternativa, que parece exigir a remoção do mapa de cores ( mencionado aqui ). Esse método também requer a divisão da varredura, por isso me faz pensar se a reclassificação original falhou devido ao uso da memória ou por algum outro motivo.
precisa saber é o seguinte
@torik Sem problemas - fico feliz em dar meus dois centavos. Eu acho que remover o mapa de cores não é o caminho a percorrer. Em vez disso, eu me concentraria na divisão de dados ou no processamento em segundo plano de 64 bits.
Aaron
@ Aaron, tendo em mente que você forneceu o código para realizar o ladrilho, como você criou o subconjunto raster usado para produzir os resultados mostrados na imagem? Concluí o ladrilhador SplitRaster (produzindo 100 subconjuntos de todo o conjunto de dados raster) e tentei fazer um loop entre todos eles para reclassificar. Infelizmente, a reclassificação falhou, resultando na mesma mensagem "Erro inesperado".
precisa saber é o seguinte
4

A whuber fez um comentário sobre o uso de ferramentas lógicas para expressar essa reclassificação . Depois de um pouco de pesquisa, descobri que o InList , como parte do conjunto de ferramentas de lógica lógica do Spatial Analyst, preencheu minha necessidade.

import arcpy

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
from arcpy.sa import InList

# Pixel values of interest, named according to Table 2 of
#  http://landcover.usgs.gov/pdf/anderson.pdf
DECIDUOUS_FOREST = 41
EVERGREEN_FOREST = 42
MIXED_FOREST = 43

inRaster = r'D:\AK_NLCD_2001_land_cover_3-13-08\ak_nlcd_2001_land_cover_3-13-08_se5.img'
accepted_raster_values = [DECIDUOUS_FOREST, EVERGREEN_FOREST, MIXED_FOREST]
filteredAlaska = InList(inRaster, accepted_raster_values)
filteredAlaska.save(r'C:\alaska\ak_woods')

É de longe a solução simplista que eu pude encontrar, executa o mais rápido e não requer consideração de agrupar o conjunto de dados original. Não há necessidade de considerar a RAM disponível da máquina, pois esta ferramenta lê diretamente do disco e armazena os resultados diretamente no disco.

Resultado filtrado do Alasca usando o InList

DoggoDougal
fonte
+1 Bem feito e uma ótima solução. Por curiosidade, quanto tempo o processamento levou?
Aaron
@ Aaron, o processamento de todo o Alasca leva 13 minutos e 23,4 segundos. O subconjunto de amostra , que é um dos 100 subconjuntos de tamanhos iguais criados por SplitRaster_management, leva 7,04 segundos.
precisa saber é o seguinte
Interessante, aproximadamente os mesmos tempos de processamento entre os dois métodos (isto é, assumindo que estávamos executando sistemas semelhantes).
Aaron
Eu tenho um Intel Core 2 Duo E6850 a 3 Ghz, 4 GB de RAM, executando o Windows 7. de 64 bits. Farei uma análise de tempo da sua solução em breve. Por enquanto, estou preso ao Arc 10.0, caso contrário, investigaria o processamento em segundo plano de 64 bits.
precisa saber é o seguinte
1

Eu usei o conjunto de dados mencionado na postagem original com uma versão 10.4 do arcmap. A reclassificação falha quando a varredura de saída é uma grade, porque as contagens de células reclassificadas estão excedendo o que pode ser armazenado no campo COUNT do IVA da grade. Quando o raster de saída é um fgdb, ele é executado com êxito em cerca de 11 minutos em uma máquina antiga de 4 núcleos executando o Windows 8. Os formatos de rasterização sem grade devem funcionar, pois usam valores flutuantes de precisão dupla para o campo de contagem. Espero que você tenha o mesmo comportamento nas versões 10.2 ou 10.3. Investigaremos o uso de um formato raster diferente para a saída padrão do Reclassify.

jt64
fonte