Python使用GDAL读写遥感图像及NDVI计算详解

南乔

发布时间：2026-10-13 | 4人浏览过

原创

本文介绍利用Python GDAL库进行遥感图像读写的方法，解析栅格矩阵与大地坐标结构。通过示例代码展示如何读取TIFF文件并计算归一化植被指数（NDVI），同时探讨实际应用中NDVI值异常及空值处理的注意事项，为遥感数据分析提供基础技术参考。

GDAL与遥感图像基础结构

GDAL作为处理空间数据的开源工具包，具备多种数据格式的读写能力。遥感图像本质上是带有大地坐标信息的栅格数据，其栅格模型主要由两部分构成：

首先是栅格矩阵，它由正方形或矩形的栅格点组成，每个点对应的数值即为像元值，在遥感影像中用于表征地物的属性。遥感影像分为单波段和多波段，其中波段代表了地物属性的种类，每个波段对应一种特定的地物属性。

其次是大地坐标，它提供了空间数据参考及地图投影信息。通过仿射矩阵，可以将行列坐标映射到具体的平面坐标上。

GDAL读写遥感数据代码实现

以下是使用GDAL进行遥感数据读写的具体代码示例：


from osgeo import gdal
import os

class GRID:

  #读图像文件
  def read_img(self,filename):
    dataset=gdal.Open(filename)    #打开文件

    im_width = dataset.RasterXSize  #栅格矩阵的列数
    im_height = dataset.RasterYSize  #栅格矩阵的行数

    im_geotrans = dataset.GetGeoTransform() #仿射矩阵
    im_proj = dataset.GetProjection() #地图投影信息
    im_data = dataset.ReadAsArray(0,0,im_width,im_height) #将数据写成数组，对应栅格矩阵

    del dataset 
    return im_proj,im_geotrans,im_data

  #写文件，以写成tif为例
  def write_img(self,filename,im_proj,im_geotrans,im_data):
    #gdal数据类型包括
    #gdal.GDT_Byte, 
    #gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
    #gdal.GDT_Float32, gdal.GDT_Float64

    #判断栅格数据的数据类型
    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
    else:
      im_bands, (im_height, im_width) = 1,im_data.shape 

    #创建文件
    driver = gdal.GetDriverByName("GTiff")      #数据类型必须有，因为要计算需要多大内存空间
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)       #写入仿射变换参数
    dataset.SetProjection(im_proj)          #写入投影

    if im_bands == 1:
      dataset.GetRasterBand(1).WriteArray(im_data) #写入数组数据
    else:
      for i in range(im_bands):
        dataset.GetRasterBand(i+1).WriteArray(im_data[i])

    del dataset

if __name__ == "__main__":
  os.chdir(r'D:\Python_Practice')            #切换路径到待处理图像所在文件夹
  run = GRID()
  proj,geotrans,data = run.read_img('LC81230402013164LGN00.tif')    #读数据
  print proj
  print geotrans
  print data
  print data.shape
  run.write_img('LC81230402013164LGN00_Rewrite.tif',proj,geotrans,data) #写数据

遥感图像NDVI计算与异常处理

在掌握GDAL遥感影像读写的基础上，可以进一步进行各类遥感图像的公式计算与统计分析。

以常见的归一化植被指数（NDVI）计算为例，只需在上述代码的倒数第二行插入相关计算代码即可实现：


import numpy as np
data = data.astype(np.float)
ndvi = (data[3]-data[2])/(data[3]+data[2])             #3为近红外波段；2为红波段
run.write_img('LC81230402013164LGN00_ndvi.tif',proj,geotrans,ndvi) #写为ndvi图像

上述计算得到的是理想状态下的NDVI。在实际处理中，NDVI值理论上应位于区间[-1, 1]内，但实际数据可能出现大于1或小于-1的情况，或者存在因坏点导致的空值nan，这需要进一步的配套处理措施。

以上内容涵盖了本文的主要知识点，希望对读者有所帮助，欢迎继续支持本站。