python script to process Landsat imagery

process_image.py

import os
import glob
import sys

"""
This script processes the raw Landsat imagery into an RGB image.
The stem of the filename is passed as the only command-line argument.
"""

if __name__ == "__main__":
    filename = sys.argv[1]

    for band in [4, 3, 2]:
        ifile = filename + "_B{0}.TIF".format(band)
        ofile = filename + "_{0}-projected.tif".format(band)
        os.system("gdalwarp -t_srs EPSG:3857 "
                  + ifile + ' ' + ofile)
        os.system("listgeo -tfw " + ofile)

    rgb_file = filename + "_RGB.tif"
    os.system("convert -combine "
              + filename + "_4-projected.tif "
              + filename + "_3-projected.tif "
              + filename + "_2-projected.tif "
              + rgb_file)

    rgb_corrected_file = filename + "_RGB-corrected.tif"
    os.system("convert -sigmoidal-contrast 6x10% "
              + rgb_file + ' '
              + rgb_corrected_file)

    os.system("cp " + filename + "_2-projected.tfw "
              + filename + "_RGB-corrected.tfw")

    rgb_projected_file = filename + "_RGB-corrected-projected.tif"
    os.system("gdalwarp -s_srs EPSG:3857 " + rgb_corrected_file
              + ' ' + rgb_projected_file)
    os.system("listgeo -tfw " + rgb_projected_file)