Rasterio.open() is only getting a single band from RGB image

Question

I am trying to combine several tiles into a single background image. I am using code from here to do the tile retrieval and merging. I've used it before on Aerial Imagery, and it worked great. Now I'm using it on Stamen Terrain tiles, and it's not working.

Specifically, it only reads one channel from RGBA tiles, and the merged image is grayscale with an inconsistent mean value.

I've pulled several adjacent Stamen terrain tiles into a folder (these look to be correct, full-color jpg and tif tiles), and then ran the merge_tiles method on them.

inputFiles = []
for name in glob.glob(inputTiles):
    if 'tif' in name:
        src = rasterio.open(name)
        array = src.read()
        print(array.shape)
        inputFiles.append(src)

mergedImage, out_trans = merge(inputFiles)

The array.shape for every tile is (1, 256, 256), but it should have 3 bands for RGB (they actually have 4! even though they are jpgs, but the 4th channel is all 1s). So the problem is not with merging itself, but in reading the Stamen tiles. Using the same code on my previous tileset still correctly reads, merges, and writes the merged tif image.

So does anybody know how to correctly read in these Stamen Terrain jpg tiles using rasterio so that all the channels are detected?

Answer 1

Those tiles are single band png with a colour table, not 3 band RGB / 4 band RGBA. (A=Alpha where 0=fully transparent and 255=fully opaque).

A colour table maps a single value to an RGB triplet. In the example stamen terrain tile below, a value of 0 in the single band is mapped to RGB 46,33,21 (actually RGBA 46,33,21,255 but all the A values in that colour table are 255).

gdalinfo <stamen tile url>/terrain/12/653/1581.png
Driver: PNG/Portable Network Graphics
Files: none associated
Size is 256, 256
...
Band 1 Block=256x1 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)
    0: 46,33,21,255
    1: 51,38,26,255
    2: 59,46,35,255
    3: 69,58,47,255
    4: 79,69,58,255
    5: 90,83,72,255
    ...
    255: 252,252,251,255

You are using a older version of the code. An updated version automatically expands the single band + colour table to 3 band RGB (relevant commit).

If you want to expand the colour table to RGB yourself, you can use something like:

gdal_translate -expand rgb input.tif output_rgb.tif
Input file size is 256, 256
0...10...20...30...40...50...60...70...80...90...100 - done.

gdalinfo output_rgb.tif 
Driver: GTiff/GeoTIFF
...
Band 1 Block=256x10 Type=Byte, ColorInterp=Red
Band 2 Block=256x10 Type=Byte, ColorInterp=Green
Band 3 Block=256x10 Type=Byte, ColorInterp=Blue

Answer 2

The answer by user2856 is the correct explanation and solution. Here is a kludge that solved the problem without me understanding how/why it was happening.

in the fetch_tile function, convert the tile image using image processing

url = tile_source.replace("{{x}}", str(x)).replace("{{y}}", str(y)).replace("{{z}}", str(z))
path = f'{temp_dir}/{x}_{y}_{z}.jpg'
urllib.request.urlretrieve(url, path)  

tileImage = plt.imread(path)
image = np.dstack([tileImage[:,:,2], tileImage[:,:,1], tileImage[:,:,0]])  
cv.imwrite(path, 255 * image)
return(path)

Now the tile is saved as a "normal" RGB image before opened by rasterio and geo-referenced. It seems that plt.imread does the work of translating the color table that allows this kludge to work.