Wie klassifiziere ich einen sehr großen Landbedeckungsdatensatz neu?

Betrachten Sie den NLCD2001 Land Cover-Datensatz für Alaska ( Download-Link ). Ich muss diesen Datensatz neu klassifizieren, damit nur Pixel mit den Werten 41, 42 und 43 erhalten bleiben. Alle anderen Pixelwerte sollten NoData werden (oder 0, falls erforderlich).

Dies scheint eine einfache Aufgabe zu sein, für die nur ein Aufruf des Reclassify-Tools erforderlich ist. Leider führt jeder Anruf zu einer vagen und nicht hilfreichen Fehlermeldung:

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)

Wie kann ich dieses Raster-Dataset neu klassifizieren? Ich verwende ArcCatalog 10.0, Build 4000, mit aktivierter Spatial Analyst-Erweiterung.

Das erste angehängte Skript hat Ihre AK NLCD-Daten in ca. 15 Minuten erfolgreich neu klassifiziert (i7, 12 GB RAM-Maschine). Da der ursprüngliche Datensatz fast 7 GB groß ist, können Speicherprobleme auftreten. Wenn Sie nicht das gesamte Dataset in einem Block verarbeiten können, teilen Sie es vor der Neuklassifizierung mit dem zweiten Skript auf. Ich empfehle, eine kleine Teilmenge der Daten zu übernehmen (Klicken Sie mit der rechten Maustaste auf Rasterebene in Inhaltsverzeichnis> Daten> Daten exportieren> Umfang (Datenrahmen) und testen Sie das erste Skript. Wenn Sie die Parameter für den Befehl "Neu klassifizieren" gewählt haben, fahren Sie mit der Neuklassifizierung fort Laden Sie alternativ das 64-Bit- Produkt für die Hintergrundverarbeitung für ArcGIS 10.1 SP1 herunter , das hier verfügbar ist . Viel Glück.

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

Bearbeiten : Wenn Sie Ihre Daten vor der Verarbeitung aufteilen müssen, sollte dieses Skript helfen:

Skript 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 machte einen Kommentar zur Verwendung logischer Werkzeuge, um diese Neuklassifizierung auszudrücken . Nach einigem Graben stellte ich fest, dass InList als Teil des Logical Math-Toolset von Spatial Analyst meinen Bedarf erfüllte.

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 ist bei weitem die einfachste Lösung, die ich finden konnte, führt die schnellste aus und erfordert keine Berücksichtigung der Kachelung des Originaldatensatzes. Der verfügbare Arbeitsspeicher des Computers muss nicht berücksichtigt werden, da dieses Tool direkt von der Festplatte liest und die Ergebnisse direkt auf der Festplatte speichert.

Ich habe den im ursprünglichen Beitrag erwähnten Datensatz mit einer 10.4 dev-Version von arcmap verwendet. Die Neuklassifizierung schlägt fehl, wenn das Ausgabe-Raster ein Raster ist, da die neu klassifizierten Zellenzahlen überlaufen, was im Feld COUNT einer Raster-Mehrwertsteuer gespeichert werden kann. Wenn das Ausgabe-Raster eine FGDB ist, wird es für mich in ungefähr 11 Minuten auf einem älteren 4-Core-Computer unter Windows 8 erfolgreich ausgeführt. Nicht-Grid-Rasterformate sollten funktionieren, da sie Gleitkommawerte mit doppelter Genauigkeit für das Zählfeld verwenden. Ich gehe davon aus, dass Sie mit den veröffentlichten Versionen 10.2 oder 10.3 dasselbe Verhalten erzielen sollten. Wir werden die Verwendung eines anderen Rasterformats für die Standardausgabe für Reclassify untersuchen.