I am trying to create a natural color (RGB) image from Landsat 8 data (11 bands). Bands 4, 3, 2, correspond to red, green, blue respectively. When I stack them, however, the images i am getting are not correct. I suspect the issue is with my astype() but I am not sure what 'type' I need.

import numpy as np
from osgeo import gdal, gdal_array

red = gdal.Open(r'/Users/myfiles/Satelite/LC08_L1TP/LC08_L1TP_039037_20130616_20170309_01_T1/LC08_L1TP_039037_20130616_20170309_01_T1_B4.TIF')
green = gdal.Open(r'/Users/myfiles/Satelite/LC08_L1TP/LC08_L1TP_039037_20130616_20170309_01_T1/LC08_L1TP_039037_20130616_20170309_01_T1_B3.TIF')
blue = gdal.Open(r'/Users/myfiles/Satelite/LC08_L1TP/LC08_L1TP_039037_20130616_20170309_01_T1/LC08_L1TP_039037_20130616_20170309_01_T1_B2.TIF')

destination = r'/Users/myfiles/Satelite/rgb.tif'

red_band = red.ReadAsArray()
green_band = green.ReadAsArray()
blue_band = blue.ReadAsArray()

rgb = np.array([red_band, green_band, blue_band])
gdal_array.SaveArray(rgb.astype('uint8'), destination, 'GTiff', red)

The above code with astype('uint8'), which I have seen recommended for others with this issue, produces the following image:


using astype(np.byte) produces a very similar image. Omitting the astype() produces this image:


which, frankly, is bizarre. I have tried creating an empty .tif and writing to it:

gtiff_driver = gdal.GetDriverByName('GTiff')

rgb = gtiff_driver.Create(destination, red.RasterXSize, red.RasterYSize, 3, gdal.GDT_Byte)



but that returns an all white raster. And using np.dstack() returns the exact same image as the first example.

I am not really sure where to go from here, or where I am going wrong. As far as I can tell I am selecting the correct RGB bands. I need the image to maintain its georeference/ would like to automate it for a number of images otherwise I would just do this in photoshop (see last image).

what i want

1 Answer 1


The issue is that array.astype() is not applying any stretching or scaling, you need to do that yourself to make the input range of 16-bit values (2^16-1 == max of 65535) fit into an 8-bit integer (2^8-1 == max of 255). You're getting integer overflow leading to that wonky image, values larger than that 8-bit max are wrapping back around from the minimum.

A quick way could be to do:

scaled = (rgb * (255 / 65535)).astype(np.uint8)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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