gdal读取多波段tif和创建新的多波段tif

from osgeo import gdal
import gdal

# 读文件
YG_dataset = gdal.Open(r'G:\数据\原始分辨率图\jiuzhaigoujingxi.tif')  # 打开文件
print(YG_dataset)
YG_width = YG_dataset.RasterXSize  # 栅格矩阵的列数
print(YG_width)
YG_height = YG_dataset.RasterYSize  # 栅格矩阵的行数
print(YG_height)
YG_geotrans = YG_dataset.GetGeoTransform()  # 仿射矩阵
print(YG_geotrans)
YG_proj = YG_dataset.GetProjection()  # 地图投影信息
print(YG_proj)
YG_data = YG_dataset.ReadAsArray(0, 0, YG_width, YG_height)  # 将数据写成数组，对应栅格矩阵
print(YG_data)
print(YG_data.shape)
b1 = YG_data[0,0:YG_height,0:YG_width]#获取第1波段
print(b1)
print(b1.shape)
b2 = YG_data[1,0:YG_height,0:YG_width]#获取第2波段
b3 = YG_data[2,0:YG_height,0:YG_width]#获取第3波段
b4 = YG_data[3,0:YG_height,0:YG_width]#获取第4波段


# 创建tif文件
driver = gdal.GetDriverByName("GTiff")
# 这里的5，5就是创建一个5x5大小的tif，后面的5是波段数
# dataset_out = driver.Create(path, cols, rows, bands, dtype)
New_YG_dataset = driver.Create(r'F:\YG\result\综合夜间灯光强度.tif', 5, 5, 5, gdal.GDT_Float32)
New_YG_dataset.SetGeoTransform(YG_geotrans)
New_YG_dataset.SetProjection(YG_proj)
print(b1[0:5, 0:5] * 1.0)
band1 = New_YG_dataset.GetRasterBand(1)
band1.WriteArray(b1[0:5, 0:5] * 1.0 + 10) # 原图是整数，乘个0.1，与浮点数可以相加
band2 = New_YG_dataset.GetRasterBand(2)
band2.WriteArray(b2[0:5, 0:5] * 1.0 + 20)
band3 = New_YG_dataset.GetRasterBand(3)
band3.WriteArray(b3[0:5, 0:5] * 1.0 + 30)
band4 = New_YG_dataset.GetRasterBand(4)
band4.WriteArray(b3[0:5, 0:5] * 1.0 + 40)
band5 = New_YG_dataset.GetRasterBand(5)
band5.WriteArray(b3[0:5, 0:5] * 1.0 + 50)
print(New_YG_dataset)
New_YG_data = New_YG_dataset.ReadAsArray(0, 0, 3, 3)  # 将数据写成数组，对应栅格矩阵
print(New_YG_data)
print(New_YG_data.shape) # (5, 3, 3)