[Gdal-dev] Landsat
Norman Vine
nhv at cape.com
Mon Mar 1 18:55:58 EST 2004
Clay, Bruce writes:
>
> Are there any functions in the GDAL package for working with Landsat
> images, in particular merging multiple layers into a single (RGB) images
> and selecting a specific subset of the image.
Probably easiest to do with Python
see attached simplistic script for merging bands to RGB and
gdal_merge.py, gdalwarp or gdal_translate for clipping ops
or extend the script todo your clipping too :-)
note attached script is realy just to demonstrate 'how' this
could be done and not really designed for production use
where you would probably want todo color correction also
HTH
Norman
-------------- next part --------------
#! /usr/bin/python
"""
Make a RGB LandSat Image with color correction
using bands 3 2 1 as R G B
Expects one argument
basename of dataset
ie without the trailing '_nnX.tif'
all 3 band files must exist in the same directory
Needs the Numeric
http://www.pfdubois.com/numpy/
the GDAL extension
http://www.remotesensing.org/gdal
prebuilt extensions included in OpenEV
http://openev.sourceforge.net/
and the Python Imaging extensions <PIL>
http://www.pythonware.com/products/pil/
Author: Norman Vine - nhv at cape.com
"""
import os,sys
from gdalnumeric import *
from PIL import Image,ImageOps
from gzip import GzipFile
# change this to suit
#DEFAULT_FILE = "p012r30_5t900908"
DEFAULT_FILE = "p004r48_4t920625"
#DEFAULT_FILE = "p012r30_5t900908"
def red_lut():
return make_lut(0.6024, -1.5, 7.0/256.0)
def green_lut():
return make_lut(1.1749, -2.8, 5.0/256.0)
def blue_lut():
return make_lut(0.8059, -1.2, 3.0/256.0)
def ColorCorrect(band, lut, clip=None):
image = ArrayToImage(band)
# transform image with our lut
image = image.point(lut)
# stretch contrast by ignoring the black background
# and clipping a bit of our histogram's tails
# original LandSat Images apparently have [0 0 0] for background
image = ImageOps.autocontrast(image,cutoff=clip,ignore=[0,254,255])
return ImageToArray(image)
"""
return ImageToArray(ImageOps.autocontrast(ArrayToImage(band).point(lut),cutoff=clip,ignore=[0]))
"""
def makeLandsatRGB(base_name, LUT=None):
""" """
print "\nProcessing %s\n"%(base_name)
if LUT is None:
print "creating LUTs\n"
LUT = [red_lut(),green_lut(),blue_lut()]
ds = getBands(base_name)
merged_ds = getMerged(base_name, ds)
BANDS = ["RED","GREEN","BLUE"]
print "\nprocessing"
for i in range(1,4,1):
print"\t band %s"%(BANDS[i-1])
print "\t\treading"
band = BandReadAsArray(ds[i-1].GetRasterBand(1),0,0,None,None)
print "\t\ttransforming"
band = ColorCorrect(band, LUT[i-1], clip=0.01)
print"\t\twriting"
BandWriteArray( merged_ds.GetRasterBand(i), band )
for i in range(1,4,1):
os.remove("%s_nn%d.tif"%(base_name,i))
print "\n=== %s.tif written ===\n"%(base_name)
def splitLandsatRGB(base_name,data_path="./"):
""" """
data_path=os.path.join(data_path,base_name)
print "\nProcessing %s\n"%(data_path)
merged_ds = gdal.Open("%s.tif"%(data_path)) #(base_name))
print "\tRed: %s_nn3.tif"%(base_name)
red_ds = gdal.GetDriverByName("GTiff").Create(
"%s_nn3.tif"%(base_name), merged_ds.RasterXSize, merged_ds.RasterYSize, 1 )
CopyDatasetInfo(merged_ds, red_ds )
BandWriteArray( red_ds.GetRasterBand(1), BandReadAsArray(merged_ds.GetRasterBand(1),0,0,None,None) )
print "\tGreen: %s_nn2.tif"%(base_name)
green_ds = gdal.GetDriverByName("GTiff").Create(
"%s_nn2.tif"%(base_name), merged_ds.RasterXSize, merged_ds.RasterYSize, 1 )
CopyDatasetInfo(merged_ds, green_ds )
BandWriteArray( green_ds.GetRasterBand(1), BandReadAsArray(merged_ds.GetRasterBand(2),0,0,None,None) )
print "\tBlue: %s_nn1.tif"%(base_name)
blue_ds = gdal.GetDriverByName("GTiff").Create(
"%s_nn1.tif"%(base_name), merged_ds.RasterXSize, merged_ds.RasterYSize, 1 )
CopyDatasetInfo(merged_ds, blue_ds )
BandWriteArray( blue_ds.GetRasterBand(1), BandReadAsArray(merged_ds.GetRasterBand(3),0,0,None,None) )
def unzip(fname):
out_file = open(fname,"wb")
zip_file = GzipFile(fname+".gz","r")
out_file.write(zip_file.read(size=-1))
def ArrayToImage(a):
i=Image.fromstring('L',(a.shape[1],a.shape[0]),
(a.astype('b')).tostring())
return i
def ImageToArray(i):
a=fromstring(i.tostring(),'b')
a.shape=i.im.size[1], i.im.size[0]
return a
def getBands(base_name):
print "unzipping"
ds = []
for i in range(3,0,-1):
getBand(ds,base_name,i)
return ds
def getBand(ds,base_name,i):
print"\t band %d"%(i)
name = "%s_nn%d.tif"%(base_name,i)
unzip(name)
ds.append(gdal.Open(name))
def getMerged(base_name, ds):
merged_ds = gdal.GetDriverByName("GTiff").Create(
"%s.tif"%(base_name), ds[0].RasterXSize, ds[0].RasterYSize, 3 )
CopyDatasetInfo(ds[0],merged_ds )
return merged_ds
def stretch(band):
mmin = minimum.reduce(ravel(band))
mmax = maximum.reduce(ravel(band))
band = (band - mmin) / (mmax - mmin)
return band
def make_lut(gain,bias,exposure):
lut = arange(0,256,1,'f')
ilut = arange(0,256,1,'b')
lut = clip(256-256*exp(-exposure*(gain*lut+bias)),0,255)
for i in range(256):
ilut[i] = lut[i]
# print ilut
return ilut
if __name__ == "__main__":
# if len(sys.argv) > 1:
# fname = sys.argv[1]
# else:
# fname = DEFAULT_FILE
# names = ["p012r30_5t900908","p011r31_5t901003","p012r31_5t870916",
# "p007r52_5t19861231","p007r53_5t19861231","p013r55_4t880206" ]
# "p004r48_4t920625","p004r47_5t910207"]
# makeLandsatRGB("p011r31_5t901003")
# names = ["p012r43_5t19840603"]
# for fname in names:
# makeLandsatRGB(dirname,fname)
# splitLandsatRGB(fname)
# os.system("../PrepData/Prepland3 %s_nn"%(fname))
# hires_names = ["1","2"]
# hires_names = ["stcroix_e","stcroix_w"]
# hires_names = ["stcroix_w"]
# hires_names = ["hadleys","quisset","juniper"]
data_path = "f:/data/processed"
for fname in hires_names:
makeLandsatRGB(fname)
splitLandsatRGB(fname,data_path)
os.system("../PrepData/Prepland3 %s_nn"%(fname))
More information about the Gdal-dev
mailing list