Comment reclasser un très grand ensemble de données sur la couverture terrestre?

Considérez l'ensemble de données sur la couverture terrestre NLCD2001 pour l'Alaska ( lien de téléchargement ). J'ai besoin de reclasser cet ensemble de données afin que seuls les pixels de valeur 41, 42 et 43 soient conservés; toutes les autres valeurs de pixels doivent devenir NoData (ou 0, si nécessaire).

Cela semble être une tâche simple, ne nécessitant qu'un seul appel à l'outil Reclassifier. Malheureusement, chaque appel entraîne un message d'erreur vague et inutile:

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)

Comment reclasser ce jeu de données raster? J'utilise ArcCatalog 10.0, Build 4000, avec l'extension Spatial Analyst activée.

Le premier script joint a reclassé avec succès vos données AK NLCD en environ 15 minutes (i7, machine de 12 Go de RAM). Étant donné que l'ensemble de données d'origine mesure près de 7 Go, vous pouvez rencontrer des problèmes de mémoire. Si vous ne pouvez pas traiter l'intégralité de l'ensemble de données en un seul bloc, essayez de le fractionner avec le deuxième script avant la reclassification. Ma recommandation est de prendre un petit sous-ensemble des données (Couche raster clic droit dans Table des matières> Données> Exporter les données> Étendue (trame de données) et tester le premier script. Une fois que vous avez composé les paramètres de la commande reclassifier, puis passer à la reclassification l'ensemble des données ou le fractionner. Sinon, essayez de télécharger le produit de géotraitement en arrière-plan 64 bits pour ArcGIS 10.1 SP1, disponible ici . Bonne chance.

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

Modifier : si vous devez diviser vos données avant le traitement, ce script devrait vous aider:

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)

whuber a fait un commentaire concernant l'utilisation d'outils logiques pour exprimer ce reclassement . Après avoir creusé un peu, j'ai trouvé InList , dans le cadre du jeu d'outils Logical Math de Spatial Analyst, qui répondait à mes besoins.

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

C'est de loin la solution la plus simple que j'ai pu trouver, s'exécute le plus rapidement et ne nécessite aucune considération de la mosaïque de l'ensemble de données d'origine. Il n'est pas nécessaire de considérer la RAM disponible de la machine, car cet outil lira directement à partir du disque et stockera les résultats directement sur le disque.

J'ai utilisé l'ensemble de données mentionné dans la publication d'origine avec une version 10,4 dev de arcmap. La reclassification échoue lorsque le raster en sortie est une grille, car le nombre de cellules reclassées déborde de ce qui peut être stocké dans le champ COUNT d'une grille TVA. Lorsque le raster en sortie est un fgdb, il s'exécute avec succès en environ 11 minutes sur une ancienne machine à 4 cœurs exécutant Windows 8. Les formats raster non-grille devraient fonctionner car ils utilisent des valeurs flottantes double précision pour le champ de comptage. Je pense que vous devriez obtenir le même comportement avec les versions 10.2 ou 10.3. Nous étudierons l'utilisation d'un format raster différent pour la sortie par défaut de Reclassification.