#载入相关库
import os
from osgeo import gdal
import matplotlib.pyplot as plt
import numpy as np

#读取数据集
os.chdir(r'D:\开发python\遥感示例数据')
band1_fn='LC08_L1TP_124039_20180517_20180517_01_RT_B3.TIF'#绿光波段
band3_fn='LC08_L1TP_124039_20180517_20180517_01_RT_B4.TIF'#红光波段
band4_fn='LC08_L1TP_124039_20180517_20180517_01_RT_B5.TIF'#近红外波段

in_ds=gdal.Open(band1_fn)
#获取dataset第一个波段（仅有一个）
in_band=in_ds.GetRasterBand(1)
gtiff_driver=gdal.GetDriverByName('GTiff')
out_ds=gtiff_driver.Create('nat_color.tiff',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_ds有三个波段，获取第三个波段，然后将in_data的数据输入第三个波段
out_band=out_ds.GetRasterBand(3)
out_band.WriteArray(in_data)
#更简略的写法
in_ds=gdal.Open(band3_fn)
out_band=out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())
#更更精简的写法
out_ds.GetRasterBand(1).WriteArray(gdal.Open(band4_fn).ReadAsArray())


out_ds.FlushCache()
for i in range(1,4):
    out_ds.GetRasterBand(i).ComputeStatistics(False)

out_ds.BuildOverviews('average',[2,4,8,16,32])

#python无法显示255之外的图像，也就是数据格式为unit16的图像，显示出来的画会是黑白图像
#img=np.dstack((out_ds.GetRasterBand(1).ReadAsArray(), out_ds.GetRasterBand(2).ReadAsArray(), out_ds.GetRasterBand(3).ReadAsArray()))
#img=out_ds.GetRasterBand(1).ReadAsArray()

#释放内存
del out_ds