Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
#!/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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.