点击上方蓝字关注我们

“ 基本操作第二辑+简单波段合成。”

一.基本操作

通过gdal提供的GetDriverByName方法构建一个driver对象。其中file_format是想要保存的数据对象。(关于gdal支持的栅格驱动对象，可以通过https://gdal.org/drivers/raster/index.html进行查询，常见的shp,hdr,tif等都包括在内)

driver.Create(filename,xsize,ysize,bands,data_type)'''filename 创建的数据集的路径xsize 新数据集的列数ysize 新数据集的行数bands 新数据集的波段数data_type 将要存储在新数据集中的数据类型'''

通过driver对象的Create方法构建一个新的数据集，通常需要填入的参数是，文件名，数据的xsize和ysize,波段通道数目，和gdal数据类型

(ps: gdal的数据类型大概有以下几种：

GDT_Byte = 1(C++中对应unsigned char)    GDT_CFloat32 = 10    GDT_CFloat64 = 11    GDT_CInt16 = 8    GDT_CInt32 = 9    GDT_Float32 = 6(C++中对应float)    GDT_Float64 = 7 (C++中对应double)    GDT_Int16 = 3(C++中对应 short 或 short int)    GDT_Int32 = 5(C++中对应int 或 long 或 long int)    GDT_TypeCount = 12    GDT_UInt16 = 2(C++中对应 unsigned short)    GDT_UInt32 = 4(C++中对应unsigned long)    GDT_Unknown = 0未知数据类型    GFT_Integer = 0    GFT_Real = 1    GFT_String = 2

)

二.简单波段合成

波段合成使用Landsat8 OLI影像，Landsat 8 是美国陆地卫星计划(Landsat)的第八颗卫星，于2013年2月11号在加利福尼亚范登堡空军基地由Atlas-V火箭搭载发射成功，最初称为“陆地卫星数据连续性任务”(Landsat Data Continuity Mission，LDCM)。Landsat 8上携带陆地成像仪(Operational Land Imager ，OLI)和热红外传感器(Thermal Infrared Sensor，TIRS)。

Landsat8 各波段介绍

选取条带号202，行列号60地区的2018年6月29日获 取的图像。通过读取Landsat8 OLI的432波段信息并将这三个波段写入同一个文件中，得到432波段合成的真彩色图像。

对于GDAL数据集，有一下的一个基本结构。

GDAL数据集的基本结构。每个数据集包含一个或多个波段，又依次包含像素数据

import osfrom osgeo import gdalos.chdir(r'F:\作业\遥感地学分析\第一次作业')  # 切换当前工作路径band1_fn = 'LC08_L1TP_201024_20180629_20180716_01_T1_B4.TIF'band2_fn = 'LC08_L1TP_201024_20180629_20180716_01_T1_B3.TIF'band3_fn = 'LC08_L1TP_201024_20180629_20180716_01_T1_B2.TIF'in_ds = gdal.Open(band1_fn)in_band = in_ds.GetRasterBand(1)gtiff_driver = gdal.GetDriverByName('GTiff') # 使用gtiff的驱动创建tif格式的文件out_ds = gtiff_driver.Create('test.tiff',                              in_band.XSize, in_band.YSize, 3, in_band.DataType) # 创建一个新的3通道数据集对象，像素值大小由原本波段提供out_ds.SetProjection(in_ds.GetProjection()) # 获取需要导入报端数据的投影信息out_ds.SetGeoTransform(in_ds.GetGeoTransform()) # 获取需要导入波段数据的geotransform信息for i, z in zip([band1_fn,band2_fn,band3_fn], [1, 2, 3]): # 打开不同波段文件，并保存波段信息     in_ds = gdal.Open(i)     in_band = in_ds.GetRasterBand(1)     in_data = in_band.ReadAsArray()     out_band = out_ds.GetRasterBand(z)     out_band.WriteArray(in_data)     print('Band%s is done' % z)out_ds.FlushCache() # 清除缓存for i in range(1, 4): # 配置金字塔     out_ds.GetRasterBand(i).ComputeStatistics(False)out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])

这里做的事主要是导入gdal模块。然后设置当前的路径，并且制定Landsat波段对应的文件。通过给gdal.Open函数传递文件名，然后打开geotiff所包含的第一个波段，获取了数据集中第一个波段的句柄(注意GetRasterBand波段数从1开始)。

接着是传递一个新的数据集将像素数据复制并放入。因此需要使用到上面提到的创建驱动对象的函数。之后从输入的数据集中得到投影并复制到新数据集，接着对geotransform做同样的操作。

创建数据集之后，进行添加像素值。通过ReadAsArray得到波段像素值，然后通过指定out_ds需要放入的波段位置，使用WriteArray将像素值复制到新的波段。

使用envi打开图像得到如下结果：

可以看到合成效果基本是真彩色的样子(这里其实可以用python可视化tiff的，但是拉伸函数过长摸了以后再写，over)