Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
    
Extract by Attributes appears to also do what I need, but unfortunately results in another "Unexpected Error". –  torik Jan 3 '13 at 16:10
    
Tried another dataset maybe? Two processes failing on the same dataset makes ya wonder... –  Chad Cooper Jan 3 '13 at 16:42
2  
Ordinarily, reclassify should be a last resort, because it is so general in scope that it likely uses methods that are less efficient than can be obtained when the reclassification is easy to express arithmetically or logically. In the present case, the criterion for reclassification is so simple you ought first to try it with Con or even straight arithmetical operations (because they are fast). For instance, "grid" * ("grid" >= 41) * ("grid" <= 43) ought to do it. RAM shouldn't be an issue--Spatial Analyst automatically windows its raster I/O and these are local operations. –  whuber Jan 4 '13 at 19:18
1  
Inlist is a nice solution (+1). I was able to use con and monitored RAM usage during the operation. It never exceeded 180 MB, which is barely greater than the RAM used just to launch ArcMap. The tiling in ArcGIS is automatic--you don't even get to control it (unless you are programming to the C/Fortran interface). It appears that RAM limitations are of little concern. –  whuber Jan 7 '13 at 18:46
1  
@whuber, con worked for me as well, with the condition "Value" >= 41 AND "Value" <= 43. I would have gone with this solution, but I'm not sure if additional raster values are going to be of interest in the future. Obviously I could add an OR into the where clause, but then it starts to become more complicated. InList seems the most straight-forward solution in regards to readability and maintainability. –  torik Jan 7 '13 at 20:42
show 5 more comments

2 Answers

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)
share|improve this answer
3  
From a performance standpoint, it would be interesting to try an alternative approach using arcpy.RasterToNumPyArray() and do the reclass in numpy. You would likely want to split the raster up into tiles anyway for memory purposes, but I know that with GDAL, re-classing numpy arrays is very fast. –  DavidF Jan 3 '13 at 20:30
    
@DavidF Agreed, there would likely be significant improvement in performance. –  Aaron Jan 3 '13 at 20:50
    
Thank you for the tips, Aaron. I will give it a run as soon as I finish another workaround, which seems to require the removal of the color map (referenced here‌​). This method requires splitting the raster as well, so it makes me wonder if Reclassify original failed due to memory-usage or some other reason. –  torik Jan 3 '13 at 22:13
    
@torik No problem--I'm happy to give my two cents. I think removing the color map is not the way to go. Rather, I would focus on splitting data or 64 bit background processing. –  Aaron Jan 4 '13 at 5:04
    
@Aaron, bearing in mind that you provided code to accomplish the tiling, how did you create the subset raster the you used to produce the pictured results? I have completed the SplitRaster tiling (producing 100 subsets of the entire raster dataset), and attempted to loop through them all to reclassify. Reclassification failed, unfortunately, resulting in the same "Unexpected Error" message. –  torik Jan 4 '13 at 15:24
show 4 more comments
up vote 3 down vote accepted

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.

Filtered Alaska result using InList

share|improve this answer
    
+1 Well done and a great solution. Out of curiousity, how long did the processing take? –  Aaron Jan 7 '13 at 18:21
    
@Aaron, processing all of Alaska takes 13 minutes and 23.4 seconds. The sample subset, which is one of the 100 equally-sized subsets created by SplitRaster_management, takes 7.04 seconds. –  torik Jan 7 '13 at 20:17
    
Interesting, roughly the same processing times between the two methods (i.e. assuming we were running similar systems). –  Aaron Jan 7 '13 at 20:30
    
I have an Intel Core 2 Duo E6850 @ 3 Ghz, 4GB of RAM, running 64-bit Windows 7. I'll do a timing analysis of your solution shortly. I'm stuck with Arc 10.0 for the time being, else I would investigate 64-bit background processing. –  torik Jan 7 '13 at 20:45
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.