¿Cómo reclasificar un conjunto de datos de cobertura terrestre muy grande?

Considere el conjunto de datos NLCD2001 Land Cover para Alaska ( enlace de descarga ). Necesito reclasificar este conjunto de datos para que solo se conserven los píxeles de valor 41, 42 y 43; todos los demás valores de píxeles deben convertirse en NoData (o 0, si es necesario).

Esto parece una tarea simple, que solo requiere una llamada a la herramienta Reclasificar. Desafortunadamente, cada llamada da como resultado un mensaje de error vago 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)

¿Cómo puedo reclasificar este dataset ráster? Estoy usando ArcCatalog 10.0, Build 4000, con la extensión Spatial Analyst habilitada.

El primer script adjunto reclasificó con éxito sus datos AK NLCD en aproximadamente 15 minutos (i7, máquina de 12 GB de RAM). Dado que el conjunto de datos original es de casi 7 GB, es posible que tenga problemas de memoria. Si no puede procesar todo el conjunto de datos en un fragmento, intente dividirlo con el segundo script antes de la reclasificación. Mi recomendación es tomar un pequeño subconjunto de los datos (haga clic con el botón derecho en la capa ráster en TOC> Datos> Exportar datos> Extensión (Marco de datos) y pruebe el primer script. Una vez que marque los parámetros para el comando de reclasificación, luego avance hacia la reclasificación todo el conjunto de datos o dividiéndolo. Alternativamente, intente descargar el producto de geoprocesamiento en segundo plano de 64 bits para ArcGIS 10.1 SP1, disponible aquí . Mucha suerte.

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 : si necesita dividir sus datos antes del procesamiento, este script debería ayudarlo:

Guión 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)

whuber hizo un comentario sobre el uso de herramientas lógicas para expresar esta reclasificación . Después de un poco de investigación, encontré que InList , como parte del conjunto de herramientas Logical Math de Spatial Analyst, llenaba mi necesidad.

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')

Es, con mucho, la solución simplista que pude encontrar, se ejecuta más rápido y no requiere considerar el mosaico del conjunto de datos original. No es necesario tener en cuenta la RAM disponible de la máquina, ya que esta herramienta leerá directamente desde el disco y almacenará los resultados directamente en el disco.

He usado el conjunto de datos mencionado en la publicación original con una versión de desarrollo 10.4 de arcmap. La reclasificación falla cuando el ráster de salida es una cuadrícula, porque los recuentos de celdas reclasificadas están desbordando lo que se puede almacenar en el campo COUNT de un IVA de cuadrícula. Cuando el ráster de salida es un fgdb, se ejecuta con éxito para mí en aproximadamente 11 minutos en una máquina de 4 núcleos más antigua que ejecuta Windows 8. Los formatos de ráster sin cuadrícula deberían funcionar ya que usan valores flotantes de precisión doble para el campo de recuento. Espero que deba tener el mismo comportamiento con las versiones 10.2 o 10.3 lanzadas. Investigaremos usando un formato ráster diferente para la salida predeterminada para Reclassify.