gdal处理landsat8合成真彩色影像tif
landsat8 影像的4，3，2波段分别对应于可见光的红光、绿光和蓝光波段，本文使用开源GIS 库jupyter python gdal 合成真彩色影像。
下载landsat8数据
在地理空间数据云上下载北京地区的landsat8影像数据
code
import os
from osgeo import gdal
import numpy as np
band1_fn='data/LC81230322021250LGN00/LC08_L1TP_123032_20210907_20210907_01_RT_B4.TIF'#red band
band2_fn='data/LC81230322021250LGN00/LC08_L1TP_123032_20210907_20210907_01_RT_B3.TIF'#green band
band3_fn='data/LC81230322021250LGN00/LC08_L1TP_123032_20210907_20210907_01_RT_B2.TIF'#blue band
#read first tif
in_ds=gdal.Open(band1_fn)
in_band=in_ds.GetRasterBand(1)
#create output tif dataset
gtiff_driver=gdal.GetDriverByName('GTiff')
out_ds=gtiff_driver.Create('data/LC81230322021250LGN00/true_color.tiff',in_band.XSize,in_band.YSize,3,in_band.DataType)
#set project and transform information
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())
#read data as array
in_data=in_band.ReadAsArray()
#read first tif and write to the first band of out_ds
out_band=out_ds.GetRasterBand(1)
out_band.WriteArray(in_data)
#read second tif and write to the second band of out_ds
in_ds=gdal.Open(band2_fn)
out_band=out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())
#read third tif and write to the third band of out_ds
out_ds.GetRasterBand(3).WriteArray(gdal.Open(band4_fn).ReadAsArray())
#save
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 memory
del out_ds
result
在QGIS中可视化：