python-gdal-处理Landsat8影像

之前处理图像的时候都是用ENVI软件进行的，一个个用鼠标点，这次学习了python-arcgis处理图像，真的是太方便了，相融合哪个波段，改几个字符就可以了。
下面展示一些代码。

#用landsat8 数据的三个波段合成一张图像融合的是5、6、4波段，主要表示陆地和水。
from osgeo import gdal
import os
os.chdir(r'F:\烨叶\work\study\ENVI\LC81160272018129LGN00\LC81160272018129LGN00')
band4_fn='LC08_L1TP_116027_20180509_20180517_01_T1_B4.tif'
band5_fn='LC08_L1TP_116027_20180509_20180517_01_T1_B5.tif'
band6_fn='LC08_L1TP_116027_20180509_20180517_01_T1_B6.tif'

in_ds=gdal.Open(band4_fn)
in_band=in_ds.GetRasterBand(1)  #使用GetRasterBand()函数时，波段数一直是从1开始的，不是从0开始的
giff_driver=gdal.GetDriverByName('GTiff')        #用GTIFF的驱动对象
#下面一句是用了Create的函数，driver.Create('fielname,xsize,ysize,[bands],[datatype],[options])
out_ds=giff_driver.Create('nat_color1.tif',in_band.XSize,in_band.YSize,3,in_band.DataType)
#了解使用的投影类型，并获得投影
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())

in_data=in_band.ReadAsArray()
out_band=out_ds.GetRasterBand(3)
out_band.WriteArray(in_data)

in_ds=gdal.Open(band6_fn)
out_band=out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())

out_ds.GetRasterBand(1).WriteArray(gdal.Open(band5_fn).ReadAsArray())

#FlushCache函数遍历波段，计算每一个波段的统计数值，False告诉他你需要的是实际的统计数据，而不是估计值
out_ds.FlushCache()
for i in range(1,4):
    out_ds.GetRasterBand(i).ComputeStatistics(False)

#建立数据集的概视图层
out_ds.BuildOverviews('average',[2,4,8,16,32])
del out_ds