How to reclassify a very large land cover dataset?
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)
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.