How to reclassify a very large land cover dataset?

Question

Consider the NLCD2001 Land Cover dataset for Alaska (download link). I need to reclassify this dataset so that only pixels of value 41, 42, and 43 are preserved; all other pixel values should become NoData (or 0, if necessary).

This seems like a simple task, only requiring one call to the Reclassify tool. Unfortunately, every call results in a vague and unhelpful error message:

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)

How can I go about reclassifying this raster dataset? I am using ArcCatalog 10.0, Build 4000, with the Spatial Analyst extension enabled.

Answer 1

The first attached script successfully reclassified your AK NLCD data in about 15 minutes (i7, 12GB RAM machine). Since the original dataset is almost 7GB you may be encountering memory issues. If you cannot process the entire dataset in one chunk, try splitting it up with the second script prior to reclassification. My recommendation is to take a small subset of the data (Right click raster layer in TOC > Data > Export Data > Extent (Data Frame) and test the first script. Once you dial in the parameters for the reclassify command, then move toward reclassifying the entire dataset or splitting it up. Alternatively, try downloading the 64 bit Background Geoprocessing product for ArcGIS 10.1 SP1, available here. Best of luck.

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

Edit: If you need to split your data up prior to processing, this script should help:

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)

Answer 2

whuber made a comment regarding the usage of logical tools to express this reclassification. After a little digging, I found InList, as part of the Logical Math toolset of Spatial Analyst, filled my need.

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

It is by far the simplist solution I could find, executes the fastest, and requires no consideration of tiling the original dataset. There is no need to consider the available RAM of the machine, as this tool will read straight from disk and store the results right back on disk.

Answer 3

I've used the dataset mentioned in the original posting with a 10.4 dev version of arcmap. Reclassify fails when the output raster is a grid, because the reclassified cell counts are overflowing what can be stored in the COUNT field of a grid VAT. When the output raster is an fgdb, it executes successfully for me in about 11 minutes on an older 4 core machine running Windows 8. Non-grid raster formats should work since they use double precision floating values for the count field. I expect you should get the same behavior with 10.2 or 10.3 released versions. We'll investigate using a different raster format for default output for Reclassify.