基于python的gdal读取遥感影像

基于python的gdal读写遥感影像

  • 1. gdal介绍
  • 2. 代码详解
  • 2.1 读取数据
  • 2.2 写入影像
  • 3. 完整案例
  • 1. gdal介绍

    GDAL(Geospatial Data Abstraction Library)主要用来读取地理空间数据,现在的GDAL包并不是单独的GDAL,而是集成了GDAL和OGR的。OGR用于处理矢量数据。因此,GDAL既可以用来处理栅格也可以处理矢量文件。

    2. 代码详解

    下面代码主要为读写tif的代码,图像其实就是一个二维矩阵,其他格式的遥感影像都可以参考。包括nc,hdf等等。

    2.1 读取数据

    # 调包
    from osgeo import gdal
    import numpy as np
    
    1. 打开影像:gdal.Open()为打开图像的操作,gdal.GetDriverByName()为注册哪种格式的数据。通常影像为Geotiff格式,也可以不写参数,默认注册所有格式。
    # 打开图像并创建空间
     data = gdal.Open(filename)
     driver = gdal.GetDriverByName('GTiff')
    
    1. 获取数据的基本信息:
    # 数据集的基本信息
        print('Raster Driver : {d}\n'.format(d=driver.ShortName))
        print('影像的波段数: ', data.RasterCount)
        img_width, img_height = data.RasterXSize, data.RasterYSize
        print('影像的列,行数: {r}rows * {c}colums'.format(r=img_width, c=img_height))
        print('栅格数据的空间参考:{}'.format(data.GetGeoTransform()))  # 栅格数据的6参数
        print('投影信息:{}\n'.format(data.GetProjection()))  # 栅格数据的投影
    

    输出如下:
    关于数据的基本信息,其中数据的空间参考投影信息最为重要。在写影像时是必不可少的参数。GetGeoTransform()为获取数据的仿射变换矩阵,6个参数,记录了影像的左上角坐标,旋转,xy方向上的分辨率;GetProjection()为获取影响的投影信息,即影像的地理参考,一般都为WGS84。
    若数据格式为NC或者HDF,仿射变换矩阵需要计算。

    Raster Driver : GTiff
    影像的波段数:  6
    影像的列行数: 1143rows * 721colums
    栅格数据的空间参考:(429255.0, 30.0, 0.0, 4429725.0, 0.0, -30.0)
    投影信息:PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]
    
    1. 将数据转换为数组形式,方便进行后期处理。
      注意:python中的索引从0开始到n-1,而波段读取的索引从1到n。
    # 读取第一个波段并将其转为array
     band_1 = data.GetRasterBand(1)
     band_1 = band_1.ReadAsArray(0, 0, img_width, img_height)  # 行,列
    

    2.2 写入影像

    写一个遥感影像分为两部分:完成这张图像,包括图像的行列号,像素点的数据类型。影像在地球上的位置,即空间信息。包括仿射变换矩阵和地理位置。

    # ①存储参数设置
    driver = gdal.GetDriverByName("GTiff")
    etype = gdal.GDT_Float32
    # ②orginal_file为原始影像,获取原始影像的空间参考和投影信息。将其写入新的影像
    proj = gdal.Open(orginal_file).GetProjection()
    transform = gdal.Open(orginal_file).GetGeoTransform()
    
    # 创建存储影像基本信息:
    # result_file结果影像
    # arr为需要存储的影像矩阵(行*列),写如的格式为先列后行。band表示波段数量,etype是数据类型
    ds = driver.Create(result_file, arr.shape[1], arr.shape[0], band, etype)        
    # eg.写单波段图像,并且将空间参考和投影信息写入。
    ds.GetRasterBand(1).WriteArray(arr)
    ds.SetGeoTransform(transform)
    ds.SetProjection(proj)
    ds.FlushCache()
    del ds
    

    3. 完整案例

    两个函数,分别为文件的读写。此处随意选择一张tif格式的数据尝试,对需要更改的参数进行了详细的说明。

    #  导入所需要的包
    from osgeo import gdal
    import numpy as np
    
    #  读取遥感影像,获取影像的基本信息
    def read_tif(filename, b=1):
    	'''
    		filename为文件名称,b为波段数量;读取遥感影像并将其每个band转化为一列变量。
    		b=1是默认参数,如果读取单波段影像,不需要输入;当输入影像为多波段时,需要将b改为影像的波段数。
    	'''
    	# 打开图像并创建空间
        data = gdal.Open(filename)
        driver = gdal.GetDriverByName('GTiff')
        # driver = data.GetDriver()
    
        # 数据集的基本信息
        print('Raster Driver : {d}\n'.format(d=driver.ShortName))
        print('影像的波段数: ', data.RasterCount)
        img_width, img_height = data.RasterXSize, data.RasterYSize
        print('影像的列,行数: {r}rows * {c}colums'.format(r=img_width, c=img_height))
        print('栅格数据的空间参考:{}'.format(data.GetGeoTransform()))  # 栅格数据的6参数
        print('投影信息:{}\n'.format(data.GetProjection()))  # 栅格数据的投影
    
        # 读取影像的元数据,获取各个波段的信息
        print(data.GetMetadata())
        band = None
        for i in range(b):
            band_1 = data.GetRasterBand(i + 1)
            band_1 = band_1.ReadAsArray(0, 0, img_width, img_height).flatten()  # 行,列
            if band is None:
                band = band_1
            else:
                band = np.vstack((band, band_1))
        return band.T
    
     def arr2raster(self, orginal_file, result_file, arr, band=1):
         '''
             orginal_file为原始影像;目的为获取原始影像的地理参考。
             result_file为要保存的文件的名称,arr为转存的矩阵,proj和transform分别为投影信息和仿射变换矩阵
             band=1,表示默认情况生成单波段图像,若要多波段图像,在输入的时候更改band为要创建的波段数
         '''
         driver = gdal.GetDriverByName("GTiff")
         etype = gdal.GDT_Float32
         proj = gdal.Open(orginal_file).GetProjection()
         transform = gdal.Open(orginal_file).GetGeoTransform()
    
         # no_data_value = ''
         if band == 1:
             ds = driver.Create(result_file, arr.shape[1], arr.shape[0], band, etype)        # 行,列
             # 写单波段图像
             ds.GetRasterBand(1).WriteArray(arr)
         else:
             ds = driver.Create(result_file, arr.shape[2], arr.shape[1], band, etype)
             # 写多波段图像
             for i in range(band):
                 ds.GetRasterBand(i + 1).WriteArray(arr[i, :, :])
         ds.SetGeoTransform(transform)
         ds.SetProjection(proj)
         ds.FlushCache()
         del ds
    
    # 案例演示:读取影像
    # 使用的是landsat5/TM影像,查看演示效果
    filename=r'C:\Users\lenovo\Desktop\test.tif'
    datasets = read_tif(filename, b=6)
    
    # 案例演示:写入影像
    # 实际中,会对影像进行操作,将结果再次转为tif。例如:将分类结果转为tif影像
    # arr为需要转存的矩阵,格式是行*列。此处将原始影像的band1转为tif,行列号721 * 1143
    arr = datasets[:,0].reshape(721, 1143)
    result_file = r'C:\Users\lenovo\Desktop\result.tif'
    arr2raster(filename, result_file, arr, band=1)
    

    原始的landsat5/TM影像,真彩色合成

    写入的单波段影像

    物联沃分享整理
    物联沃-IOTWORD物联网 » 基于python的gdal读取遥感影像

    发表评论