Friday, 10 May 2019

arcpy - How to reclassify a very large land cover dataset?


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



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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...