Getting RGB values of each pixel in a LANDSAT 8 GeoTIFF with Python 3

Question

I'm working with a set of LANDSAT 8 Collection 1 Tier 1 GeoTIFFS pulled from Google Earth Engine, and I need to find the RGB values of every pixel in this GeoTIFF.

I've tried using the code from this answer to a similar question:

from osgeo import gdal
gdal.UseExceptions()

ds = gdal.Open(fname)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
ct = band.GetColorTable()

# index value to RGB (ignore A)
i2rgb = [ct.GetColorEntry(i)[:3] for i in range(ct.GetCount())] 

But I'm getting the following error:

Traceback (most recent call last):
  File "get_color.py", line 12, in <module>
    i2rgb = [ct.GetColorEntry(i)[:3] for i in range(ct.GetCount())]
AttributeError: 'NoneType' object has no attribute 'GetCount'

Ideally, I'd like to have a matrix of (R, G, B) tuples or lists, but I'm open to other formats.

Here's my gdalinfo dump:

Driver: GTiff/GeoTIFF
Files: LC08_014030_20170220.tif
Size is 749, 1570
Coordinate System is:
PROJCS["WGS 84 / UTM zone 18N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                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",-75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32618"]]
Origin = (604200.000000000000000,4855080.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  604200.000, 4855080.000) ( 73d42'13.61"W, 43d50'29.23"N)
Lower Left  (  604200.000, 4807980.000) ( 73d42'46.30"W, 43d25' 2.88"N)
Upper Right (  626670.000, 4855080.000) ( 73d25'27.68"W, 43d50'16.59"N)
Lower Right (  626670.000, 4807980.000) ( 73d26' 7.41"W, 43d24'50.42"N)
Center      (  615435.000, 4831530.000) ( 73d34' 8.82"W, 43d37'40.10"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Description = B1
Band 2 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B2
Band 3 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B3
Band 4 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B4
Band 5 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B5
Band 6 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B6
Band 7 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B7
Band 8 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B8
Band 9 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B9
Band 10 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B10
Band 11 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B11
Band 12 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = BQA

Answer 1

Your Landsat raster doesn't have a colour table, it's raw data. If you actually want the values from bands 4,3,2 which correspond to the red, green and blue visible wavelengths, you can just read as a numpy NDArray.

import numpy as np
from osgeo import gdal

ds = gdal.Open(your_geotiff)
rgb = np.stack((ds.GetRasterBand(b).ReadAsArray() for b in (4,3,2)))

Answer 2

I have not seen the GetColorTable() function in GDAL, but you could:

extract each band like:

    band1 = ds.GetRasterBand(1).ReadAsArray()
    band2 = ds.GetRasterBand(2).ReadAsArray()
    band3 = ds.GetRasterBand(3).ReadAsArray()

and stack the arrays with numpy:

    import numpy as np
    rgb_array = np.dstack([band1,band2,band3])

this produces a n x m x 3 ndarray