Create RGB images from landsat scenes. Defualts to band pattern for true color images from landsat8, Requires https://github.com/gina-alaska/dans-gdal-scripts

composite_landsat.py

#!/usr/bin/env python
from osgeo import gdal
import os
import sys
import glob
from optparse import OptionParser

def convertTo8Bit(infile, outfile):
    os.system('gdal_contrast_stretch -ndv 0 -percentile-range 0.02 0.98 %s %s' %  (infile, outfile))


def compositeDir(source_dir, dest_file, bands, pan_band, lum_bands):
    downsampled_files = []
    rgb_string = ''
    for band in bands.split(','):
        try:
            band_filename = glob.glob("%s/*B%s.TIF" % (source_dir, band))[0]
            (band_filename, downsampled) = downSampleIfNeeded(band_filename)
            if downsampled:
                downsampled_files.append(band_filename)
            rgb_string += ' -rgb ' + band_filename
        except Exception, e:
            print e

    lum_string = ''
    for l in lum_bands.split(','):
        (band, percent) = l.split(':')
        band_filename = glob.glob("%s/*B%s.TIF" % (source_dir, band))[0]
        (band_filename, downsampled) = downSampleIfNeeded(band_filename)
        if downsampled:
            downsampled_files.append(band_filename)
        lum_string += ' -lum %s %s ' % (band_filename, percent)

    pan_band_filename = glob.glob("%s/*B%s.TIF" % (source_dir, pan_band))[0]
    (pan_band_filename, downsampled) = downSampleIfNeeded(pan_band_filename)
    if downsampled:
        downsampled_files.append(pan_band_filename)

    merge_command = 'gdal_landsat_pansharp %s %s -pan %s -ndv 0 -o %s' % (rgb_string, lum_string, pan_band_filename, outfile)
    print(merge_command)
    os.system(merge_command)
    for to_remove in downsampled_files:
        os.remove(to_remove)


def downSampleIfNeeded(filename):
    ds = gdal.Open(filename, gdal.GA_ReadOnly)
    if ds is None:
        print 'Could not open ' + filename
        return (None, False)
    
    hDriver = ds.GetDriver();
    #print( "Driver: %s/%s" % ( hDriver.ShortName, hDriver.LongName ))
    metadata = papszMetadata = ds.GetMetadata_List()
    
    structure_metadata = ds.GetMetadata_List("IMAGE_STRUCTURE")
    if structure_metadata is None:
        return (None, False)

    band = ds.GetRasterBand(1)
    dataType = gdal.GetDataTypeName(band.DataType)
    
    if dataType is None:
        return (None, False)

    if dataType != "Byte":
        print "Image is " + str(nbits) + " bits, downsampling"
        outfile = filename + ".8.tiff"
        if not os.path.exists(outfile):
            print "downsampled file already exists"
            convertTo8Bit(filename, outfile)
        return (outfile, True)
        
    return (filename, False)


if __name__=='__main__':
    usage = "usage: %prog -b 4,3,2 imagedir"
    parser = OptionParser(usage=usage,
        description='''Convert a multi band landsat scene to a single RGB image,
defaults are set to create a true color composite from landsat 8 images''')
    parser.add_option("-b", "--bands", action="store", type="string", dest="bands",
     default='4,3,2')
    parser.add_option("-p", "--panchromatic", action="store", type="string", dest="pan_band",
     default='8')
    parser.add_option("-l", "--lum", action="store", type="string", dest="lum_bands",
     default='2:0.25,3:0.23,4:0.52')
    parser.add_option("-o", "--outfile", action="store", type="string", dest="outfile")

    (options, args) = parser.parse_args()

    gdal.AllRegister()

    for dir in args:
        outfile = options.outfile
        if outfile is None:
            outfile = os.path.join(dir, 'composite-%s.tiff' % options.bands)
        compositeDir(dir, outfile, options.bands, options.pan_band, options.lum_bands)