利用python批量从img影像中提取所需的波段

￡神魔★判官ぃ 2022-11-18 04:40 394阅读 0赞

**利用python批量从img影像中提取所需的波段**  
现在我知道landsat数据下载，得到的波段是一张张单个的影像，但是MODIS影像下载出来的波段是合在一起的，如果处理的数据量很大，那么手动将其波段信息提取出来很麻烦，所以接下来的代码，可以将img影像中所需要的波段提取出来，废话不多说，看代码。

**1.注意事项：**

注：这个代码要用高版本的python，本人这个是3.7的，而且如果没有gdal numpy osgeo库需要自行配置，如果没有这几个库，代码无法运行成功。
    from osgeo import gdal  
    import numpy as np
    import os

**2.原代码如下：**  
我所需要的是红绿蓝的波段

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxODI4ODE4_size_16_color_FFFFFF_t_70]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxODI4ODE4_size_16_color_FFFFFF_t_70 1]

# -*- coding: utf-8 -*-
    import glob
    
    from osgeo import gdal
    import numpy as np
    import os
    
    '''
    说明：
        提取img数据中的红、绿、近红外波段数据
            Rrs_555	绿波段（62）
            Rrs_645	红波段（63）
            Rrs_859	近红外波段（64）
            注：下划线后的数字为对应波段的中心波长
    '''
    
    def writeTiff(im_data, im_width, im_height, im_bands, im_geotrans, im_proj, path):  # 定义函数writeTiff
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32
    
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        elif len(im_data.shape) == 2:
            im_data = np.array([im_data])
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape
        # 创建文件
        driver = gdal.GetDriverByName("GTiff")
        dataset = driver.Create(path, im_width, im_height, im_bands, datatype)
        if (dataset != None):
            dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
            dataset.SetProjection(im_proj)  # 写入投影
        for i in range(im_bands):
            dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
        del dataset
    
    
    if __name__ == '__main__':
        gdal.AllRegister()
    
        path_mother = r'C:\Users\Administrator\Desktop\1' # 输入路径，文件夹
        path_out = r'C:\Users\Administrator\Desktop\3'  # 输出路径，文件夹
    
        # 获取该目录下所有文件，存入列表中
        os.chdir(path_mother)
        for file in glob.glob('*.img'):
            img_file_name = path_mother + os.sep + file
            dataset = gdal.Open(img_file_name)
    
            filename = file.split('.')[0]
            adfGeoTransform = dataset.GetGeoTransform()
    
            band_1 = dataset.GetRasterBand(62)   # 读取数据集的属性，选择待提取波段.62,63,64代表的是波段其在img中的序列号（绿波段在img影像中位于第62个波段）
            band_2 = dataset.GetRasterBand(63)
            band_3 = dataset.GetRasterBand(64)
    
            im_width = dataset.RasterXSize  # 栅格矩阵的列数
            im_height = dataset.RasterYSize  # 栅格矩阵的行数
    
            Rrs_1 = band_1.ReadAsArray(0, 0, im_width, im_height)  # 目标波段
            Rrs_2 = band_2.ReadAsArray(0, 0, im_width, im_height)
            Rrs_3 = band_3.ReadAsArray(0, 0, im_width, im_height)
    
            im_geotrans = dataset.GetGeoTransform()  # 获取仿射矩阵信息
            im_proj = dataset.GetProjection()  # 获取投影信息
    
            im_bands = 1  # band_469.RasterCount  # 输出文件的波段数
            B1 = Rrs_1 * 1.0  # / 10000
            B2 = Rrs_2 * 1.0  # / 10000
            B3 = Rrs_3 * 1.0  # / 10000
    
            outpath_1 = path_out + os.sep + filename + '_Rrs_555.Green' + '.tif'
            writeTiff(B1, im_width, im_height, im_bands, im_geotrans, im_proj, outpath_1)
    
            outpath_2 = path_out + os.sep + filename + '_Rrs_645.Red' + '.tif'
            writeTiff(B2, im_width, im_height, im_bands, im_geotrans, im_proj, outpath_2)
    
            outpath_3 = path_out + os.sep + filename + '_Rrs_859.NIR' + '.tif'
            writeTiff(B3, im_width, im_height, im_bands, im_geotrans, im_proj, outpath_3)
            print('this is ok')
            print('------All is ok------')
       
     希望以上代码能够解决大家的问题！
     长路漫漫....
     唯有坚持....

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxODI4ODE4_size_16_color_FFFFFF_t_70]: /images/20221022/83f355f3fb2d4683add3977d833890b0.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQxODI4ODE4_size_16_color_FFFFFF_t_70 1]: /images/20221022/c82abb8f4c5d4c909abe888b18f0370e.png