1

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
  • Thanks for the fast reply @Luke. I've added the gdalinfo dump to the original post. I'm a bit new to GIS processing, so if there's anything else that would be helpful please let me know and I'll get it in my post as soon as possible. – Pierogi11 Sep 13 '18 at 0:27
0

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)))
0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.