Jak przeklasyfikować bardzo duży zestaw danych dotyczących pokrycia terenu?

Zastanów się nad zestawem danych dotyczących pokrycia terenu NLCD2001 dla Alaski ( link do pobrania ). Muszę przeklasyfikować ten zestaw danych, aby zachowane były tylko piksele o wartości 41, 42 i 43; wszystkie pozostałe wartości w pikselach powinny stać się NoData (lub 0, jeśli to konieczne).

To wydaje się prostym zadaniem, wymagającym tylko jednego połączenia z narzędziem Przeklasyfikowania. Niestety każde połączenie powoduje niejasny i nieprzydatny komunikat o błędzie:

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)

Jak mogę zmienić klasyfikację tego zestawu danych rastrowych? Korzystam z ArcCatalog 10.0, Build 4000, z włączonym rozszerzeniem Spatial Analyst.

Pierwszy dołączony skrypt pomyślnie przeklasyfikował dane AK NLCD w ciągu około 15 minut (i7, 12 GB pamięci RAM). Ponieważ oryginalny zestaw danych ma prawie 7 GB, możesz mieć problemy z pamięcią. Jeśli nie możesz przetworzyć całego zestawu danych w jednym fragmencie, spróbuj podzielić go na drugi skrypt przed przeklasyfikowaniem. Polecam wziąć niewielki podzbiór danych (kliknij prawym przyciskiem myszy warstwę rastrową w Spisie treści> Dane> Eksportuj dane> Rozszerz (ramka danych) i przetestuj pierwszy skrypt. Po wybraniu parametrów polecenia przeklasyfikowania, przejdź do przeklasyfikowania cały zestaw danych lub jego podział. Alternatywnie spróbuj pobrać 64-bitowy produkt Geoprocessing w tle dla ArcGIS 10.1 SP1, dostępny tutaj . Powodzenia.

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

Edycja : jeśli musisz podzielić swoje dane przed przetwarzaniem, ten skrypt powinien pomóc:

Skrypt 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 skomentował użycie narzędzi logicznych do wyrażenia tej przeklasyfikowania . Po krótkim kopaniu znalazłem InList , jako część zestawu narzędzi Logical Math programu Spatial Analyst, który spełnił moją potrzebę.

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

Jest to zdecydowanie najprostsze rozwiązanie, jakie udało mi się znaleźć, działa najszybciej i nie wymaga uwzględnienia oryginalnego zestawu danych. Nie ma potrzeby rozważania dostępnej pamięci RAM maszyny, ponieważ to narzędzie będzie czytać bezpośrednio z dysku i zapisywać wyniki z powrotem na dysku.

Użyłem zestawu danych wspomnianego w oryginalnym poście z wersją arcmap w wersji 10.4 dev. Przeklasyfikowanie kończy się niepowodzeniem, gdy wyjściowym rastrem jest siatka, ponieważ liczba przeklasyfikowanych komórek jest przepełniona, co można zapisać w polu COUNT podatku VAT od siatki. Kiedy wyjściowy raster jest fgdb, wykonuje się dla mnie z powodzeniem za około 11 minut na starszym 4-rdzeniowym komputerze z systemem Windows 8. Nie-gridowe formaty rastrowe powinny działać, ponieważ używają zmiennoprzecinkowych wartości podwójnej precyzji dla pola zliczania. Oczekuję, że powinieneś uzyskać takie samo zachowanie w wersjach wydanych 10.2 lub 10.3. Zbadamy użycie innego formatu rastrowego jako domyślnego wyjścia dla Przeklasyfikowania.