Using NASA Above landcover imagery in place with ArcGIS

Landsat-derived Annual Dominant Land Cover Across ABoVE Core Domain is composite tiled classified imagery dataset covering the years 1984-2014. The format is elegant and clever but hard to use. Each tile file is a 31-band GeoTIFF, with each band providing one year of data for the period (band_1 = 1984, …band_31 = 2014) and classified into 15 values (1-15 in integer).

Neither ArcGIS or Qgis can use the data very well out of the box as they only allow colour maps and/or attributes for single channel 8 bit rasters. We could explode the bands into distinct files, but for my area of interest that will yield 175*31=5,425 files. A management headache to say the least, and probably file-system performance limit as well not to mention the increase in storage needs.

My intuition says both Arc and Qgis (via gdal) have the necessary features to be able to effectively use these imagery in place via raster mosaic-dataset with function chains (arc) and virtual raster tables vrt with pixel functions (gdal, and arc to lesser degree). I have so far though not managed to figure out the necessary magic.

The end goal is to enable people to load a given year from the data source and have it displayed in the standard colour map. Ideally it also yield the corresponding code when interrogated in ArcMap/Pro or Qgis. It should also be straightforward to use an arbitraty year as input to a geoprocessing tool. All with the least amount of fuss and hairpulling. 😉

I would prefer a single solution that can serve both ArcGIS and Qgis but as requested I am splitting this into a separate question for each platform. This is the ArcGIS posting. The qgis posting is here.

For ArcGIS the simplest process for a single tile I’ve devised so far is below. Unfortunately it doesn’t work for raster mosaic datasets.

• Use Catalog to explore to the rasters, expand to see the bands
  • Add a band to a map. Band_5 would be 1989.
  • Change Symbology to "Unique Values"
    • Refer to a table for the cell (pixel) value classifications and/or
    • Import canned colour map from a layer file

I’m going to self-answer with the best I’ve come up with so far for mosaic solutions, but they’re just partial steps along the way and I’m not happy with them. Please feel free to improve or better.

import arcpy
coordsys="PROJCS['Albers_Conic_Equal_Area',GEOGCS['GCS_GRS_1980_IUGG_1980',DATUM['D_unknown',SPHEROID['GRS80',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['false_easting',0.0],PARAMETER['false_northing',0.0],PARAMETER['central_meridian',-96.0],PARAMETER['standard_parallel_1',50.0],PARAMETER['standard_parallel_2',70.0],PARAMETER['latitude_of_origin',40.0],UNIT['Meter',1.0]]"

band_definitions = ";".join([f"y{x} 1 15" for x in range(1984,2015)])
    # creates "y1984 1 15; y1985 1 15; ..." which is: 
    #   {band name} {min value} {max value};

with arcpy.EnvManager(outputCoordinateSystem=coordsys):
    arcpy.management.CreateMosaicDataset(r"Z:ExternalNasanasa_imagery.gdb", 
        "src_above",
        coordsys,
        31, # number of bands
        "8_BIT_UNSIGNED", 
        "CUSTOM", 
        band_definitions)

arcpy.management.AddRastersToMosaicDataset("src_above", 
    "Raster Dataset", 
    r"Z:ExternalNasaAbove", # source folder path 
    "UPDATE_CELL_SIZES", 
    "UPDATE_BOUNDARY", 
    "NO_OVERVIEWS", 
    None, 
    0, 1500, 
    None, 
    r"*.tif", # pattern match these files only
    "NO_SUBFOLDERS", 
    "OVERWRITE_DUPLICATES", 
    "NO_PYRAMIDS", 
    "NO_STATISTICS", 
    "NO_THUMBNAILS", 
    '', 
    "NO_FORCE_SPATIAL_REFERENCE", 
    "NO_STATISTICS",
     None, 
     "NO_PIXEL_CACHE")