Export GeoTIFF with color table from web service using gdal_translate

Question

The gdal_translate program has an option to define an XML file with the necessary information to download an image using various web map server specifications. I am trying to download a file that should have a specific color table palette with 256 colors associated with it and one band of 8-bit (byte) values. But when I run gdal_translate on an XML file and download the file, the file only appears to have 3-band RGB values with no color table (when I look at the image attributes using gdalinfo).

For example, I have the following XML file I'm using to download an image:

<GDAL_WMS>
     <Service name="AGS">
         <ServerUrl>https://raster.nationalmap.gov/arcgis/rest/services/LandCover/USGS_EROS_LandCover_NLCD/MapServer</ServerUrl>
         <ImageFormat>tiff</ImageFormat>
         <Layers>show:6</Layers>
         <BBoxOrder>xyXY</BBoxOrder>
         <SRS>26917</SRS>
     </Service>
     <DataWindow>
         <UpperLeftX>720586.0306923158900</UpperLeftX>
         <UpperLeftY>3962270.5642631231000</UpperLeftY>
         <LowerRightX>723891.5076759960500</LowerRightX>
         <LowerRightY>3960158.2215603264000</LowerRightY>
         <SizeX>164</SizeX>
         <SizeY>105</SizeY>
     </DataWindow>
     <UnsafeSSL>true</UnsafeSSL>
</GDAL_WMS>

I save this file as "nlcd.xml". Then I run the following command:

gdal_translate -of GTiff "F:\temp\nlcd.xml" "f:\temp\r.nlcd.tif"

The file is downloaded and I am able to view the image using various methods. However, when I run gdalinfo on the image, there are 3 bands with no color table defined:

gdalinfo f:\temp\r.nlcd.tif

Driver: GTiff/GeoTIFF
Files: f:\temp\r.nlcd.tif
Size is 164, 105
Coordinate System is:
PROJCS["NAD83 / UTM zone 17N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-81],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","26917"]]
Origin = (720586.030692315890000,3962270.564263123100000)
Pixel Size = (20.155347461464373,-20.117549550444597)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  720586.031, 3962270.564) ( 78d33'34.26"W, 35d46'47.19"N)
Lower Left  (  720586.031, 3960158.222) ( 78d33'36.36"W, 35d45'38.69"N)
Upper Right (  723891.508, 3962270.564) ( 78d31'22.72"W, 35d46'44.50"N)
Lower Right (  723891.508, 3960158.222) ( 78d31'24.85"W, 35d45'36.00"N)
Center      (  722238.769, 3961214.393) ( 78d32'29.55"W, 35d46'11.60"N)
Band 1 Block=164x16 Type=Byte, ColorInterp=Red
Band 2 Block=164x16 Type=Byte, ColorInterp=Green
Band 3 Block=164x16 Type=Byte, ColorInterp=Blue

I expected to see one band with a color table defined-something like the following on a file I previously downloaded using another method:

gdalinfo "F:\temp\nlcd2006.tif"

Driver: GTiff/GeoTIFF
Files: F:\temp\nlcd2006.tif
Size is 1296, 1598
Coordinate System is:
PROJCS["WGS 84 / UTM zone 10N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235604902,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-123],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32610"]]
Origin = (464430.000000000000000,4363320.000000000000000)
Pixel Size = (30.000000000000000,-30.012515644555695)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=PACKBITS
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  464430.000, 4363320.000) (123d24'47.61"W, 39d25' 7.21"N)
Lower Left  (  464430.000, 4315360.000) (123d24'38.52"W, 38d59'11.44"N)
Upper Right (  503310.000, 4363320.000) (122d57'41.57"W, 39d25' 9.83"N)
Lower Right (  503310.000, 4315360.000) (122d57'42.41"W, 38d59'14.02"N)
Center      (  483870.000, 4339340.000) (123d11'12.52"W, 39d12'11.43"N)
Band 1 Block=1296x1 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)
    0: 0,0,0,255
    1: 0,250,0,255
    2: 0,0,0,255
    3: 0,0,0,255
    4: 0,0,0,255
    5: 0,0,0,255
    6: 0,0,0,255
    7: 0,0,0,255
    8: 0,0,0,255
    9: 0,0,0,255
   10: 0,0,0,255
   11: 71,107,161,255
   12: 209,222,250,255
   13: 0,0,0,255
   14: 0,0,0,255
   15: 0,0,0,255
   16: 0,0,0,255
   17: 0,0,0,255
   18: 0,0,0,255
   19: 0,0,0,255
   20: 0,0,0,255
   21: 222,201,201,255
   22: 217,148,130,255
   23: 237,0,0,255
   24: 171,0,0,255
   25: 0,0,0,255
   26: 0,0,0,255
   27: 0,0,0,255
   28: 0,0,0,255
   29: 0,0,0,255
   30: 0,0,0,255
   31: 179,173,163,255
   32: 250,250,250,255
...

My question is: Is there a way to export a GeoTIFF with a single band and a color palette from a web service using gdal_translate instead of exporting the RGB values in 3 separate bands?

Answer 1

You can now accomplish this with rasterio. From the example at: https://rasterio.readthedocs.io/en/latest/topics/color.html

import rasterio

with rasterio.Env():

    with rasterio.open('tests/data/shade.tif') as src:
        shade = src.read(1)
        meta = src.meta

    with rasterio.open('/tmp/colormap.tif', 'w', **meta) as dst:
        dst.write(shade, indexes=1)
        dst.write_colormap(
            1, {
                0: (255, 0, 0, 255),
              255: (0, 0, 255, 255) })
        cmap = dst.colormap(1)

where you would replace the example's dictionary -- which only has entries for indexes 0 and 255 -- with your own colormap dictionary.

I have found that you can supply a dictionary with 3-member tuples (R, G, B) and rasterio will give the transparency ("A") channel a value of 255.

Answer 2

I figured out a workable solution to this problem after thinking about it for a couple of months. Maybe there's an easier solution, but I have not found one through my searching and asking around. The solution I decided on was to go ahead and download the RGB image using gdal_translate and the XML file as I described in the original question. Then I create a simple text file that lists all the possible RGB values in a downloaded land use file and the index values and names I want to use in the new raster file. For example, I created the following file for NLCD 2006 data after getting the list of possible land use ID's/colors from the full original dataset that I downloaded here:

ID R G B "Name"
---------------
0 0 0 0 "Unclassified"
11 71 107 160 "Open Water"
12 209 221 249 "Perennial Snow/Ice"
21 221 201 201 "Developed, Open Space"
22 216 147 130 "Developed, Low Intensity"
23 237 0 0 "Developed, Medium Intensity"
24 170 0 0 "Developed, High Intensity"
31 178 173 163 "Barren Land"
41 104 170 99 "Deciduous Forest"
42 28 99 48 "Evergreen Forest"
43 181 201 142 "Mixed Forest"
52 204 186 124 "Shrub/Scrub"
71 226 226 193 "Herbaceuous"
81 219 216 61 "Hay/Pasture"
82 170 112 40 "Cultivated Crops"
90 186 216 234 "Woody Wetlands"
95 112 163 186 "Emergent Herbaceuous Wetlands"

Then I wrote some C++ code to read the downloaded RGB GeoTIFF as a temporary file and convert this file to a file containing the desired color index values along with the color table palette. I also write a raster attribute table since I may use/need this at some point in the future. So here's my basic code that uses GDAL classes to read the temporary RGB file and convert it to an index file with a color palette if you're interested:

//------------------------------------------------------------------------------
/// \brief   Converts an RGB image to an image with a color table.
//------------------------------------------------------------------------------
void iConvertRgbToPaletteImage(const std::string &a_tempFile,
  const std::string &a_fileName, const std::string &a_rasterAttFilename)
{
  GDALAllRegister();
  auto inRaster = (GDALDataset *)GDALOpen(a_tempFile.c_str(), GA_ReadOnly);
  if (inRaster != nullptr)
  {
    std::vector<GDALRasterBand *> rBand = { inRaster->GetRasterBand(1), inRaster->GetRasterBand(2),
      inRaster->GetRasterBand(3) };
    int nCols(0);
    int nRows(0);
    std::vector<Uchar *> scanLines;
    for (const auto &band : rBand)
    {
      nCols = band->GetXSize();
      nRows = band->GetYSize();
      auto scanLine = (Uchar *)CPLMalloc(sizeof(Uchar)*nCols*nRows);
      band->RasterIO(GF_Read, 0, 0, nCols, nRows, scanLine, nCols, nRows, GDALDataType::GDT_Byte,
        0, 0);
      scanLines.push_back(scanLine);
    }
    auto attClrs = iReadRasterAttFile(a_rasterAttFilename);
    const int size(nRows * nCols);
    std::vector<Uchar> vals(size, 0);
    for (int i = 0; i < size; ++i)
    {
      auto it = attClrs.find(Color(scanLines[0][i], scanLines[1][i], scanLines[2][i]));
      if (it != attClrs.end()) vals[i] = it->second.first;
    }
    for (const auto &scanLine : scanLines) CPLFree(scanLine);
    auto driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    if (driver == nullptr) return;
    auto outRaster = driver->Create(a_fileName.c_str(), nCols, nRows, 1, GDALDataType::GDT_Byte, nullptr);
    double transform[6] = { 0, 0, 0, 0, 0, 0 };
    auto err = inRaster->GetGeoTransform(transform);
    if (err != CPLErr::CE_None) return;
    err = outRaster->SetGeoTransform(transform);
    if (err != CPLErr::CE_None) return;
    err = outRaster->SetProjection(inRaster->GetProjectionRef());
    GDALClose((GDALDatasetH)inRaster);
    std::remove(a_tempFile.c_str());
    if (err != CPLErr::CE_None) return;
    auto band = outRaster->GetRasterBand(1);
    band->SetColorInterpretation(GDALColorInterp::GCI_PaletteIndex);
    GDALColorTable t;
    vector<string> types;
    const int maxUchar(256);
    for (int i = 0; i < maxUchar; ++i)
    {
      auto it = attClrs.begin();
      for (; it != attClrs.end(); ++it)
      {
        if (it->second.first == i) break;
      }
      GDALColorEntry e;
      if (it != attClrs.end())
      {
        e.c1 = it->first.r;
        e.c2 = it->first.g;
        e.c3 = it->first.b;
        e.c4 = maxUchar;
        types.push_back(it->second.second);
      }
      else
      {
        e.c1 = e.c2 = e.c3 = e.c4 = 0;
        types.push_back("");
      }
      t.SetColorEntry(i, &e);
    }
    band->SetColorTable(&t);
    GDALDefaultRasterAttributeTable rat;
    err = rat.InitializeFromColorTable(&t);
    const int fieldIdx = rat.GetColumnCount();
    err = rat.CreateColumn("Land Cover Class", GDALRATFieldType::GFT_String,
      GDALRATFieldUsage::GFU_Name);
    if (err == CPLErr::CE_None)
    {
      for (int i = 0; i < maxUchar; ++i)
      {
        rat.SetValue(i, fieldIdx, types[i].c_str());
      }
    }
    band->SetDefaultRAT(&rat);
    band->RasterIO(GDALRWFlag::GF_Write, 0, 0, nCols, nRows,
      (void *)&vals[0], nCols, nRows, GDALDataType::GDT_Byte, 0, 0);
    GDALClose((GDALDatasetH)outRaster);
  }
} // iConvertRgbToPaletteImage

Note that "Color" is just a color class we use with a red, green, and blue value stored as unsigned chars (Uchars). There's a function in this code called "iReadRasterAttFile" that reads the attribute file shown above. Anybody who uses this code would have to figure out the best way for them to store/read this raster attribute file, so I won't include our program-specific code for this here (it is actually very simple). The prototype for this function is as follows:

std::map<Color, std::pair<Uchar, std::string>> iReadRasterAttFile(const std::string
  &a_rasterAttFilename);

This solution has worked for my needs. Does anybody have a better solution or is there an easier way to do this? I would be interested in knowing.