在Python中使用GDAL库，相当于拥有了一个功能强大的“地理数据翻译官”，可以让你用Python代码轻松处理和转换海量的卫星影像、高程数据等地理空间信息。
下面我从安装配置、核心对象、常见操作这几个方面，为你梳理在Python中使用GDAL的要点。
💻 安装与配置
在Python中安装GDAL最推荐的方法是使用 conda，它可以自动帮你解决复杂的依赖问题，避免很多麻烦。
安装方式 推荐指数 适用场景与说明 示例命令
使用 conda ⭐⭐⭐ 最推荐。尤其适合Anaconda用户，能自动处理GDAL的底层依赖库（如PROJ、GEOS等），避免编译错误。 conda install -c conda-forge gdal
使用 pip ⭐⭐ 如果你的环境是纯Python且不想用conda，可以尝试。但在Windows上可能会遇到编译问题。 pip install gdal 或借助 gdal-installer 工具简化安装
系统包管理器 ⭐⭐ 适用于Ubuntu/Debian等Linux系统，命令简单，但可能安装的版本较旧。 sudo apt-get install python3-gdal
🔍 核心对象与功能
在Python的 osgeo 模块中，有几个核心对象是你需要熟悉的：
gdal.Open(): 这是一切操作的起点，用来打开一个栅格文件（如GeoTIFF、HDF等），返回一个Dataset对象。
Dataset对象: 代表打开的整个栅格文件。通过它可以获取栅格的宽度 (RasterXSize)、高度 (RasterYSize)、波段数 (RasterCount)、地理坐标信息 (GetGeoTransform()) 和投影信息 (GetProjection())。
Band对象: 通过 dataset.GetRasterBand(i) 获取，代表栅格数据中的一个波段（例如Landsat影像中的红波段）。你可以用 band.ReadAsArray() 将整个波段的数据读取为一个NumPy数组，这是进行科学计算和分析的关键一步。反过来，也可以用 band.WriteArray() 将NumPy数组的数据写回栅格。
ogr模块: 用于处理矢量数据（如Shapefile、GeoJSON）。ogr.Open()用于打开矢量文件。
🚀 常见应用场景与代码片段
1. 读取并查看栅格数据
这是最基础的操作，你可以快速了解一个栅格文件的基本信息，并获取其像素值数组
from osgeo import gdalimport numpy as np# 打开栅格文件dataset = gdal.Open('your_image.tif')# 获取基本信息print("栅格尺寸: {} x {} (宽 x 高)".format(dataset.RasterXSize, dataset.RasterYSize))print("波段数: ", dataset.RasterCount)print("投影: ", dataset.GetProjection())# 读取第一个波段的数据为NumPy数组band1 = dataset.GetRasterBand(1)data_array = band1.ReadAsArray().astype(np.float32)print("数据数组的形状: ", data_array.shape)# 操作完毕后，将数据集设为None以释放资源dataset = None
2. 按矢量边界裁剪栅格
这是GIS分析中非常常见的需求。通常的思路是，先获取矢量边界的地理范围，然后利用这个范围去裁剪栅格
from osgeo import gdal, ogr# 打开矢量文件，获取其空间范围shp_ds = ogr.Open('your_shapefile.shp')lyr = shp_ds.GetLayer()# ... (此处省略了根据矢量范围计算裁剪窗口的代码) ...# 打开要裁剪的栅格input_raster = gdal.Open('input_raster.tif')# 使用gdal.Warp函数进行裁剪# 关键是要设置cutlineDSName和cropToCutline参数# output_raster = gdal.Warp('output_cropped.tif', input_raster,# cutlineDSName='your_shapefile.shp',# cropToCutline=True, dstNodata=0)print("裁剪操作已完成。")# 更详细的实现请参考官方文档或相关博客[citation:1][citation:4]
3. 创建新栅格并写入数组
当你通过NumPy计算出新的结果（比如NDVI植被指数）后，需要将其保存为一个带地理参考的新GeoTIFF文件
from osgeo import gdal, osrimport numpy as np# 假设我们已经有一个计算结果数组 new_datanew_data = np.random.rand(500, 500) # 示例数据# 获取驱动driver = gdal.GetDriverByName('GTiff')# 创建一个新文件，需要指定文件名、宽、高、波段数、数据类型output_dataset = driver.Create('output_ndvi.tif', xsize=new_data.shape[1], ysize=new_data.shape[0], bands=1, eType=gdal.GDT_Float32)# 将数组写入新文件的第一个波段output_band = output_dataset.GetRasterBand(1)output_band.WriteArray(new_data)# 设置地理变换参数和投影（通常从一个已有的参考文件复制）# 假设我们有一个参考数据集 ref_dataset# output_dataset.SetGeoTransform(ref_dataset.GetGeoTransform())# output_dataset.SetProjection(ref_dataset.GetProjection())# 设置波段无数据值output_band.SetNoDataValue(-9999)# 保存并关闭文件output_dataset = None
总的来说，在Python中，GDAL是你处理遥感数据最忠实的伙伴。熟练掌握 gdal.Open、ReadAsArray、WriteArray 和 gdal.Warp 这几个核心操作，就能应对绝大部分日常工作。
