简介本资源是一份面向GIS开发人员、遥感数据处理工程师及地理信息专业学习者的Python实战指南聚焦GDAL库在Geotiff格式遥感影像读写中的核心应用。内容系统讲解如何使用osgeo.gdal完成影像元数据地理变换参数、坐标系投影提取、多波段像素数据高效读取以及基于NumPy数组生成带地理参考的Geotiff文件覆盖从数据加载、空间信息保持到结果导出的完整闭环流程适用于遥感图像预处理、地理配准、算法验证等典型场景。资源为单文件PDF文档1个33KB结构紧凑含可直接复用的read_tiff与array2raster函数实现、关键API说明及调用示例代码注释清晰适合作为开发速查手册或教学补充材料。目前已有5297人学习下载对初学者掌握GDAL基础IO操作、避免常见坐标丢失或数据类型错误具有明确指导价值。1. 用 Python GDAL 打开 GeoTIFF 遥感影像不是“读个文件”那么简单你双击一个.tif文件Windows 可能用照片查看器打开——但那只是像素GIS 软件能显示坐标、波段、投影、像元大小是因为它读出了嵌在文件头里的地理空间元数据。GDALGeospatial Data Abstraction Library正是这套“读懂遥感影像语言”的底层引擎。Python 通过osgeo.gdal模块调用它才能真正实现按地理范围裁剪、按波段组合生成真彩色图、把经纬度坐标转成图像行列号、批量重采样多景 Landsat 影像、甚至导出带坐标的 NumPy 数组用于深度学习训练。这不是cv2.imread()或PIL.open()能做到的——它们看不见 WGS84、UTM、仿射变换参数也读不出GCPs地面控制点。本文面向已装好 Python 的遥感处理初学者和 GIS 开发者不讲编译源码不绕过 conda/pip 安装坑直接从gdal.Open()第一行开始拆解如何可靠读取、验证、写入、调试 GeoTIFF尤其聚焦国产高分、Sentinel-2、Landsat 等常见遥感数据的实际参数配置。2. 安装与验证避开 Windows 下 GDAL 的 DLL 冲突和 Linux 的 pkg-config 缺失GDAL 是 C/C 库Python 绑定需二进制兼容。直接pip install gdal极易失败——它默认编译本地版本而多数用户缺gcc、proj头文件或libtiff-dev。必须用预编译轮子且版本需对齐。2.1 推荐安装路径跨平台实测有效提示不要用pip install GDAL首字母大写正确包名是小写gdal但 pip 会自动标准化关键是版本号必须与系统环境匹配。Windows推荐 conda# 创建干净环境避免与已有 Python 冲突 conda create -n gdal-env python3.9 conda activate gdal-env # 从 conda-forge 安装含完整依赖proj、geos、libtiff conda install -c conda-forge gdal3.8.4验证命令from osgeo import gdal print(gdal.__version__) # 输出 3.8.4 ds gdal.Open(test.tif) # 不报错即基础可用Ubuntu/Debian避免 apt-get 的老旧版本# 先装系统级依赖关键否则 pip 编译必跪 sudo apt update sudo apt install -y libproj-dev libgeos-dev libtiff-dev libpng-dev libjpeg-dev # 再用 pip 安装预编译 wheel指定最新稳定版 pip install --find-links https://download.osgeo.org/gdal/3.8.4/ --no-index gdal3.8.4macOSM1/M2 芯片注意架构# 使用 miniforge非 anaconda避免 x86 兼容问题 brew install miniforge conda activate base conda install -c conda-forge gdal3.8.42.2 验证 GeoTIFF 读取能力的三步检查法仅gdal.Open()成功不代表能读遥感数据。必须验证三个核心能力2.2.1 检查驱动是否识别为 GTifffrom osgeo import gdal ds gdal.Open(landsat8_b4.tif) driver_name ds.GetDriver().ShortName print(driver_name) # 必须输出 GTiff若为 MEM 或 VRT 则文件损坏或路径错误注意GetDriver().ShortName返回驱动名不是文件扩展名。有些.tif实际是 BigTIFF 或 Cloud Optimized GeoTIFFCOGGDAL 3.5 才完全支持 COG 的流式读取。2.2.2 提取关键元数据字段遥感处理刚需# 获取地理参考六参数仿射变换 geotransform ds.GetGeoTransform() print(GeoTransform:, geotransform) # 输出示例: (120.1, 0.0002695, 0.0, 30.5, 0.0, -0.0002695) # 含义: (左上角经度, 经度方向像元宽, 旋转项, 左上角纬度, 纬度方向旋转项, 纬度方向像元高) # 获取投影字符串WKT 格式 proj ds.GetProjection() print(Projection WKT length:, len(proj)) # 非空且长度 100 才可信 # 示例开头: PROJCS[WGS 84 / UTM zone 50N,GEOGCS[WGS 84,DATUM... # 获取波段数和数据类型 band_count ds.RasterCount data_type gdal.GetDataTypeName(ds.GetRasterBand(1).DataType) print(f波段数: {band_count}, 数据类型: {data_type}) # 常见: UInt16, Float322.2.3 读取第一波段统计值快速判断数据有效性band ds.GetRasterBand(1) stats band.GetStatistics(False, True) # (approx_ok, force) if stats is not None: print(fMin: {stats[0]:.2f}, Max: {stats[1]:.2f}, Mean: {stats[2]:.2f}) else: print(未计算统计值尝试 band.ComputeStatistics(False)) band.ComputeStatistics(False) stats band.GetStatistics(False, True)提示GetStatistics()返回None很常见——GeoTIFF 默认不内嵌统计值。此时必须调用ComputeStatistics(False)强制计算耗时但必要。若ComputeStatistics后仍为空说明波段数据损坏或内存不足。3. 读取与写入从单波段灰度到多波段真彩色图的完整链路读写 GeoTIFF 的核心是gdal.Dataset和gdal.Band对象。重点不是“怎么读”而是“读什么”和“怎么写对”。3.1 读取按需加载避免 OOM内存溢出遥感影像动辄 GB 级。ReadAsArray()直接全载入内存是新手最大误区。3.1.1 分块读取Tile-based Readingimport numpy as np def read_tiff_tile(filepath, xoff0, yoff0, xsizeNone, ysizeNone, band_listNone): 按偏移量和尺寸读取子区域band_list 指定波段索引1-based ds gdal.Open(filepath) if xsize is None: xsize ds.RasterXSize if ysize is None: ysize ds.RasterYSize if band_list is None: band_list list(range(1, ds.RasterCount 1)) # 预分配数组(波段数, 行, 列) data np.empty((len(band_list), ysize, xsize), dtypenp.float32) for i, bidx in enumerate(band_list): band ds.GetRasterBand(bidx) # ReadAsArray(xoff, yoff, xsize, ysize) data[i] band.ReadAsArray(xoff, yoff, xsize, ysize).astype(np.float32) return data, ds.GetGeoTransform(), ds.GetProjection() # 示例读取左上角 1000x1000 像元区域所有波段 tile_data, gt, proj read_tiff_tile(sentinel2_b04_b03_b02.tif, xoff0, yoff0, xsize1000, ysize1000) print(Tile shape:, tile_data.shape) # (3, 1000, 1000)参数说明xoff/yoff是图像左上角起始列/行号整数xsize/ysize是读取宽度/高度。ReadAsArray()内部自动处理数据类型转换但显式.astype()可控精度。3.1.2 按地理范围裁剪Ground Coordinate Croppingdef crop_by_extent(filepath, extent, epsg4326): extent: (min_lon, min_lat, max_lon, max_lat) —— WGS84 经纬度 epsg: 输入 extent 的坐标系默认 4326 ds gdal.Open(filepath) gt ds.GetGeoTransform() proj ds.GetProjection() # 将地理范围转为图像行列号 from osgeo import osr srs_in osr.SpatialReference() srs_in.ImportFromEPSG(epsg) srs_out osr.SpatialReference() srs_out.ImportFromWkt(proj) transform osr.CoordinateTransformation(srs_in, srs_out) # 转换四个角点 corners [ transform.TransformPoint(extent[0], extent[1]), # min_lon, min_lat transform.TransformPoint(extent[2], extent[1]), # max_lon, min_lat transform.TransformPoint(extent[2], extent[3]), # max_lon, max_lat transform.TransformPoint(extent[0], extent[3]), # min_lon, max_lat ] # 取并集行列范围 xs [int((c[0] - gt[0]) / gt[1]) for c in corners] ys [int((c[1] - gt[3]) / gt[5]) for c in corners] xoff, xsize min(xs), max(xs) - min(xs) 1 yoff, ysize min(ys), max(ys) - min(ys) 1 return read_tiff_tile(filepath, xoff, yoff, xsize, ysize) # 裁剪北京五环内区域WGS84 beijing_extent (116.2, 39.7, 116.5, 40.0) crop_data, _, _ crop_by_extent(gf1_pms1.tiff, beijing_extent)注意CoordinateTransformation要求输入坐标系与目标坐标系均有效。若proj为空字符串需先用ds.SetProjection()设置否则TransformPoint报错。3.2 写入生成带完整地理信息的新 GeoTIFF写入的关键是复现GeoTransform和Projection否则输出只是普通 TIFF。3.2.1 创建单波段 GeoTIFF如 NDVI 结果图def write_geotiff(filepath, array, geotransform, projection, nodataNone, dtypegdal.GDT_Float32): array: 2D numpy 数组行, 列 geotransform: 6元组同 GetGeoTransform() 输出 projection: WKT 字符串同 GetProjection() 输出 driver gdal.GetDriverByName(GTiff) # 创建数据集行、列、波段数、数据类型 ds driver.Create(filepath, array.shape[1], array.shape[0], 1, dtype) ds.SetGeoTransform(geotransform) ds.SetProjection(projection) band ds.GetRasterBand(1) if nodata is not None: band.SetNoDataValue(nodata) band.WriteArray(array) # 必须关闭数据集否则文件不写入磁盘 ds.FlushCache() del ds # 显式释放 # 示例将计算好的 NDVI 写入 ndvi_array (nir_band.astype(float) - red_band.astype(float)) / (nir_band red_band 1e-8) write_geotiff(ndvi_result.tif, ndvi_array, gt, proj, nodata-9999)3.2.2 创建多波段 GeoTIFF真彩色合成def write_rgb_geotiff(filepath, r_array, g_array, b_array, geotransform, projection): r/g/b_array: 均为 2D numpy 数组形状相同 driver gdal.GetDriverByName(GTiff) ds driver.Create(filepath, r_array.shape[1], r_array.shape[0], 3, gdal.GDT_UInt16) ds.SetGeoTransform(geotransform) ds.SetProjection(projection) # 写入 R、G、B 波段1-indexed ds.GetRasterBand(1).WriteArray(r_array) ds.GetRasterBand(2).WriteArray(g_array) ds.GetRasterBand(3).WriteArray(b_array) # 设置波段描述可选但利于后续识别 ds.GetRasterBand(1).SetDescription(Red Band) ds.GetRasterBand(2).SetDescription(Green Band) ds.GetRasterBand(3).SetDescription(Blue Band) ds.FlushCache() del ds # 合成 Sentinel-2 真彩色B04Red, B03Green, B02Blue write_rgb_geotiff(sentinel2_truecolor.tif, b04, b03, b02, gt, proj)参数说明Create()第四参数是数据类型gdal.GDT_UInt16适配 12bit 遥感数据若用Float32需确保array是 float 类型否则静默截断。4. 常见故障排查为什么读出来全是 0为什么写入后 QGIS 打不开生产环境中80% 的 GDAL 问题源于元数据缺失、类型不匹配或路径权限。4.1 读取全零或异常值的 5 个检查点检查项命令/代码期望结果问题定位1. 文件是否被其他进程锁定lsof | grep your_file.tif(Linux/macOS) 或资源监视器 (Windows)无输出文件被 ENVI/QGIS 占用GDAL 只读模式失败2. 波段数据类型是否匹配print(band.DataType)1(Byte),2(UInt16),6(Float32)若误用ReadAsArray()读UInt16到int8数组高位被截断为 03. NoData 值是否被忽略print(band.GetNoDataValue())如0,65535,-9999读取后需array[array nodata] np.nan4. 图像压缩是否不支持print(ds.GetMetadata(IMAGE_STRUCTURE))COMPRESSION: LZW或DEFLATEGDAL 3.4 支持 LZW/DEFLATE旧版需--config GDAL_PAM_ENABLED NO5. 内存映射是否启用gdal.Info(file.tif, formatjson)[metadata][IMAGE_STRUCTURE]INTERLEAVE: BAND或PIXELBAND交错存储更省内存PIXEL适合 OpenCV 处理4.1.1 快速诊断脚本保存为gdal_check.pyfrom osgeo import gdal import sys def check_tiff(filepath): ds gdal.Open(filepath) if ds is None: print(f❌ 打开失败{filepath}) return print(f✅ 文件打开成功) print(f 驱动: {ds.GetDriver().ShortName}) print(f 尺寸: {ds.RasterXSize}x{ds.RasterYSize}x{ds.RasterCount}) gt ds.GetGeoTransform() print(f GeoTransform: {gt[:4]}... (前4位)) proj ds.GetProjection() print(f 投影长度: {len(proj)} 字符) for i in range(1, min(3, ds.RasterCount 1)): b ds.GetRasterBand(i) dt gdal.GetDataTypeName(b.DataType) ndv b.GetNoDataValue() print(f 波段{i}: {dt}, NoData{ndv}) # 读取左上角 3x3 像元检查值 win b.ReadAsArray(0, 0, 3, 3) print(f 左上角值: {win.flatten()}) if __name__ __main__: check_tiff(sys.argv[1])运行python gdal_check.py landsat8_b5.tif输出清晰显示每层状态比gdalinfo更聚焦遥感处理痛点。4.2 写入后 QGIS/GIS 软件打不开的硬核修复QGIS 报错 “Invalid layer” 或 “Projection not recognized” 通常因 WKT 投影字符串不标准。4.2.1 强制标准化投影解决 EPSG 代号丢失from osgeo import osr def standardize_projection(wkt_str): 将任意 WKT 投影转为标准 EPSG 代号格式 srs osr.SpatialReference() srs.ImportFromWkt(wkt_str) # 尝试导出为 EPSG 代号如存在 epsg_code srs.AutoIdentifyEPSG() if epsg_code 0: return wkt_str # 无法识别返回原字符串 else: srs_out osr.SpatialReference() srs_out.ImportFromEPSG(int(srs.GetAuthorityCode(None))) return srs_out.ExportToWkt() # 修复后写入 standard_proj standardize_projection(proj) write_geotiff(fixed.tif, array, gt, standard_proj)4.2.2 添加地理参考当 GeoTransform 为 (0,1,0,0,0,1) 时def fix_identity_geotransform(ds, lon_min, lat_min, lon_max, lat_max): 用地理范围反推仿射变换适用于无地理参考的 TIFF xsize, ysize ds.RasterXSize, ds.RasterYSize xres (lon_max - lon_min) / xsize yres (lat_max - lat_min) / ysize # 左上角坐标 gt (lon_min, xres, 0, lat_max, 0, -yres) ds.SetGeoTransform(gt) # 设置 WGS84 投影 srs osr.SpatialReference() srs.ImportFromEPSG(4326) ds.SetProjection(srs.ExportToWkt()) return ds # 示例给一张纯 RGB 图添加 WGS84 地理参考 ds gdal.Open(rgb.jpg, gdal.GA_Update) ds fix_identity_geotransform(ds, 116.2, 39.7, 116.5, 40.0) ds.FlushCache() del ds注意gdal.GA_Update模式打开才能写入元数据。若原文件只读需先shutil.copy()到新路径再操作。5. 进阶技巧用 GDAL 处理国产高分影像的坐标偏移与波段顺序国产高分系列GF-1/2/6/7影像常含两个特殊问题1WGS84 坐标系下存在百米级偏移2波段顺序与国际标准不一致如 GF-1 PMS 的蓝绿红近红外顺序是 B2,B3,B4,B5而非 B1,B2,B3,B4。GDAL 提供原生接口修正。5.1 用 GCPs地面控制点校正坐标偏移高分影像的.xml附带 GCP 列表GDAL 可直接读取并生成 RPC 模型。5.1.1 从 XML 提取 GCPs 并写入 GeoTIFFimport xml.etree.ElementTree as ET def parse_gf_gcps(xml_path): 解析高分 XML 中的 GCPs tree ET.parse(xml_path) root tree.getroot() gcps [] for gcp in root.findall(.//GCP): line float(gcp.find(Line).text) samp float(gcp.find(Sample).text) lat float(gcp.find(Lat).text) lon float(gcp.find(Lon).text) height float(gcp.find(Height).text) if gcp.find(Height) is not None else 0 gcps.append(gdal.GCP(lon, lat, height, samp, line)) return gcps # 将 GCPs 写入 TIFF替代 GeoTransform ds gdal.Open(gf1_pms1.tif, gdal.GA_Update) gcps parse_gf_gcps(GF1_PMS1_E116.3_N39.8_20230501_L1A0000111111.xml) ds.SetGCPs(gcps, WGS84) # 生成 .aux.xml 辅助文件QGIS 识别 GCP ds.BuildOverviews(NEAREST, [2,4,8]) ds.FlushCache() del ds效果QGIS 加载后自动启用 GCP 校正消除系统性偏移。BuildOverviews生成金字塔加速大图浏览。5.2 自动识别高分波段并重排顺序高分数据常以B*命名波段但 GDAL 默认按文件内顺序读取。需解析元数据确定物理波长。5.2.1 读取波段中心波长并排序def get_band_wavelengths(ds): 从 GDAL 元数据或 XML 解析各波段中心波长nm # 方法1检查 GDAL 元数据 meta ds.GetMetadata() if BAND_WAVELENGTH in meta: return [float(x) for x in meta[BAND_WAVELENGTH].split(,)] # 方法2解析附带 XML高分标准 import glob xml_files glob.glob(ds.GetDescription().replace(.tif, *.xml)) if xml_files: tree ET.parse(xml_files[0]) wavelengths [] for band in tree.findall(.//Band): wl band.find(CenterWavelength) if wl is not None: wavelengths.append(float(wl.text)) return wavelengths # 默认波长GF-1 PMS return [450, 520, 630, 830] # B2,B3,B4,B5 def reorder_bands_to_rgb(ds, target_order[2, 1, 0]): # B4(R), B3(G), B2(B) 按波长重排波段生成真彩色数组 wavelengths get_band_wavelengths(ds) # 按波长升序排序索引蓝绿红近红外 sorted_idx sorted(range(len(wavelengths)), keylambda i: wavelengths[i]) # 映射到目标顺序取第2个红、第1个绿、第0个蓝 rgb_idx [sorted_idx[i] for i in target_order] data np.stack([ ds.GetRasterBand(i1).ReadAsArray().astype(np.float32) for i in rgb_idx ], axis0) return data # 使用 ds gdal.Open(GF1_PMS1.tif) rgb_data reorder_bands_to_rgb(ds) # 自动匹配 B4,B3,B2 write_rgb_geotiff(gf1_truecolor.tif, *rgb_data, ds.GetGeoTransform(), ds.GetProjection())逻辑说明target_order[2,1,0]表示取波长第三长红、第二长绿、第一长蓝的波段。sorted_idx是按波长从小到大排列的原始波段索引rgb_data保证输出符合人眼感知。至此你已掌握用 Python GDAL 处理 GeoTIFF 遥感影像的全链路从避坑安装、分块读取、地理裁剪、带参写入到国产数据偏移校正与波段重排。下一步可延伸至gdal.Warp投影转换、gdal.Translate格式转换、或接入rasterio做更 Pythonic 的操作——但核心原理始终是GeoTransform、Projection、RasterBand这三根支柱。本文还有配套的精品资源点击获取
简介本资源是一份面向GIS开发人员、遥感数据处理工程师及地理信息专业学习者的Python实战指南聚焦GDAL库在Geotiff格式遥感影像读写中的核心应用。内容系统讲解如何使用osgeo.gdal完成影像元数据地理变换参数、坐标系投影提取、多波段像素数据高效读取以及基于NumPy数组生成带地理参考的Geotiff文件覆盖从数据加载、空间信息保持到结果导出的完整闭环流程适用于遥感图像预处理、地理配准、算法验证等典型场景。资源为单文件PDF文档1个33KB结构紧凑含可直接复用的read_tiff与array2raster函数实现、关键API说明及调用示例代码注释清晰适合作为开发速查手册或教学补充材料。目前已有5297人学习下载对初学者掌握GDAL基础IO操作、避免常见坐标丢失或数据类型错误具有明确指导价值。1. 用 Python GDAL 打开 GeoTIFF 遥感影像不是“读个文件”那么简单你双击一个.tif文件Windows 可能用照片查看器打开——但那只是像素GIS 软件能显示坐标、波段、投影、像元大小是因为它读出了嵌在文件头里的地理空间元数据。GDALGeospatial Data Abstraction Library正是这套“读懂遥感影像语言”的底层引擎。Python 通过osgeo.gdal模块调用它才能真正实现按地理范围裁剪、按波段组合生成真彩色图、把经纬度坐标转成图像行列号、批量重采样多景 Landsat 影像、甚至导出带坐标的 NumPy 数组用于深度学习训练。这不是cv2.imread()或PIL.open()能做到的——它们看不见 WGS84、UTM、仿射变换参数也读不出GCPs地面控制点。本文面向已装好 Python 的遥感处理初学者和 GIS 开发者不讲编译源码不绕过 conda/pip 安装坑直接从gdal.Open()第一行开始拆解如何可靠读取、验证、写入、调试 GeoTIFF尤其聚焦国产高分、Sentinel-2、Landsat 等常见遥感数据的实际参数配置。2. 安装与验证避开 Windows 下 GDAL 的 DLL 冲突和 Linux 的 pkg-config 缺失GDAL 是 C/C 库Python 绑定需二进制兼容。直接pip install gdal极易失败——它默认编译本地版本而多数用户缺gcc、proj头文件或libtiff-dev。必须用预编译轮子且版本需对齐。2.1 推荐安装路径跨平台实测有效提示不要用pip install GDAL首字母大写正确包名是小写gdal但 pip 会自动标准化关键是版本号必须与系统环境匹配。Windows推荐 conda# 创建干净环境避免与已有 Python 冲突 conda create -n gdal-env python3.9 conda activate gdal-env # 从 conda-forge 安装含完整依赖proj、geos、libtiff conda install -c conda-forge gdal3.8.4验证命令from osgeo import gdal print(gdal.__version__) # 输出 3.8.4 ds gdal.Open(test.tif) # 不报错即基础可用Ubuntu/Debian避免 apt-get 的老旧版本# 先装系统级依赖关键否则 pip 编译必跪 sudo apt update sudo apt install -y libproj-dev libgeos-dev libtiff-dev libpng-dev libjpeg-dev # 再用 pip 安装预编译 wheel指定最新稳定版 pip install --find-links https://download.osgeo.org/gdal/3.8.4/ --no-index gdal3.8.4macOSM1/M2 芯片注意架构# 使用 miniforge非 anaconda避免 x86 兼容问题 brew install miniforge conda activate base conda install -c conda-forge gdal3.8.42.2 验证 GeoTIFF 读取能力的三步检查法仅gdal.Open()成功不代表能读遥感数据。必须验证三个核心能力2.2.1 检查驱动是否识别为 GTifffrom osgeo import gdal ds gdal.Open(landsat8_b4.tif) driver_name ds.GetDriver().ShortName print(driver_name) # 必须输出 GTiff若为 MEM 或 VRT 则文件损坏或路径错误注意GetDriver().ShortName返回驱动名不是文件扩展名。有些.tif实际是 BigTIFF 或 Cloud Optimized GeoTIFFCOGGDAL 3.5 才完全支持 COG 的流式读取。2.2.2 提取关键元数据字段遥感处理刚需# 获取地理参考六参数仿射变换 geotransform ds.GetGeoTransform() print(GeoTransform:, geotransform) # 输出示例: (120.1, 0.0002695, 0.0, 30.5, 0.0, -0.0002695) # 含义: (左上角经度, 经度方向像元宽, 旋转项, 左上角纬度, 纬度方向旋转项, 纬度方向像元高) # 获取投影字符串WKT 格式 proj ds.GetProjection() print(Projection WKT length:, len(proj)) # 非空且长度 100 才可信 # 示例开头: PROJCS[WGS 84 / UTM zone 50N,GEOGCS[WGS 84,DATUM... # 获取波段数和数据类型 band_count ds.RasterCount data_type gdal.GetDataTypeName(ds.GetRasterBand(1).DataType) print(f波段数: {band_count}, 数据类型: {data_type}) # 常见: UInt16, Float322.2.3 读取第一波段统计值快速判断数据有效性band ds.GetRasterBand(1) stats band.GetStatistics(False, True) # (approx_ok, force) if stats is not None: print(fMin: {stats[0]:.2f}, Max: {stats[1]:.2f}, Mean: {stats[2]:.2f}) else: print(未计算统计值尝试 band.ComputeStatistics(False)) band.ComputeStatistics(False) stats band.GetStatistics(False, True)提示GetStatistics()返回None很常见——GeoTIFF 默认不内嵌统计值。此时必须调用ComputeStatistics(False)强制计算耗时但必要。若ComputeStatistics后仍为空说明波段数据损坏或内存不足。3. 读取与写入从单波段灰度到多波段真彩色图的完整链路读写 GeoTIFF 的核心是gdal.Dataset和gdal.Band对象。重点不是“怎么读”而是“读什么”和“怎么写对”。3.1 读取按需加载避免 OOM内存溢出遥感影像动辄 GB 级。ReadAsArray()直接全载入内存是新手最大误区。3.1.1 分块读取Tile-based Readingimport numpy as np def read_tiff_tile(filepath, xoff0, yoff0, xsizeNone, ysizeNone, band_listNone): 按偏移量和尺寸读取子区域band_list 指定波段索引1-based ds gdal.Open(filepath) if xsize is None: xsize ds.RasterXSize if ysize is None: ysize ds.RasterYSize if band_list is None: band_list list(range(1, ds.RasterCount 1)) # 预分配数组(波段数, 行, 列) data np.empty((len(band_list), ysize, xsize), dtypenp.float32) for i, bidx in enumerate(band_list): band ds.GetRasterBand(bidx) # ReadAsArray(xoff, yoff, xsize, ysize) data[i] band.ReadAsArray(xoff, yoff, xsize, ysize).astype(np.float32) return data, ds.GetGeoTransform(), ds.GetProjection() # 示例读取左上角 1000x1000 像元区域所有波段 tile_data, gt, proj read_tiff_tile(sentinel2_b04_b03_b02.tif, xoff0, yoff0, xsize1000, ysize1000) print(Tile shape:, tile_data.shape) # (3, 1000, 1000)参数说明xoff/yoff是图像左上角起始列/行号整数xsize/ysize是读取宽度/高度。ReadAsArray()内部自动处理数据类型转换但显式.astype()可控精度。3.1.2 按地理范围裁剪Ground Coordinate Croppingdef crop_by_extent(filepath, extent, epsg4326): extent: (min_lon, min_lat, max_lon, max_lat) —— WGS84 经纬度 epsg: 输入 extent 的坐标系默认 4326 ds gdal.Open(filepath) gt ds.GetGeoTransform() proj ds.GetProjection() # 将地理范围转为图像行列号 from osgeo import osr srs_in osr.SpatialReference() srs_in.ImportFromEPSG(epsg) srs_out osr.SpatialReference() srs_out.ImportFromWkt(proj) transform osr.CoordinateTransformation(srs_in, srs_out) # 转换四个角点 corners [ transform.TransformPoint(extent[0], extent[1]), # min_lon, min_lat transform.TransformPoint(extent[2], extent[1]), # max_lon, min_lat transform.TransformPoint(extent[2], extent[3]), # max_lon, max_lat transform.TransformPoint(extent[0], extent[3]), # min_lon, max_lat ] # 取并集行列范围 xs [int((c[0] - gt[0]) / gt[1]) for c in corners] ys [int((c[1] - gt[3]) / gt[5]) for c in corners] xoff, xsize min(xs), max(xs) - min(xs) 1 yoff, ysize min(ys), max(ys) - min(ys) 1 return read_tiff_tile(filepath, xoff, yoff, xsize, ysize) # 裁剪北京五环内区域WGS84 beijing_extent (116.2, 39.7, 116.5, 40.0) crop_data, _, _ crop_by_extent(gf1_pms1.tiff, beijing_extent)注意CoordinateTransformation要求输入坐标系与目标坐标系均有效。若proj为空字符串需先用ds.SetProjection()设置否则TransformPoint报错。3.2 写入生成带完整地理信息的新 GeoTIFF写入的关键是复现GeoTransform和Projection否则输出只是普通 TIFF。3.2.1 创建单波段 GeoTIFF如 NDVI 结果图def write_geotiff(filepath, array, geotransform, projection, nodataNone, dtypegdal.GDT_Float32): array: 2D numpy 数组行, 列 geotransform: 6元组同 GetGeoTransform() 输出 projection: WKT 字符串同 GetProjection() 输出 driver gdal.GetDriverByName(GTiff) # 创建数据集行、列、波段数、数据类型 ds driver.Create(filepath, array.shape[1], array.shape[0], 1, dtype) ds.SetGeoTransform(geotransform) ds.SetProjection(projection) band ds.GetRasterBand(1) if nodata is not None: band.SetNoDataValue(nodata) band.WriteArray(array) # 必须关闭数据集否则文件不写入磁盘 ds.FlushCache() del ds # 显式释放 # 示例将计算好的 NDVI 写入 ndvi_array (nir_band.astype(float) - red_band.astype(float)) / (nir_band red_band 1e-8) write_geotiff(ndvi_result.tif, ndvi_array, gt, proj, nodata-9999)3.2.2 创建多波段 GeoTIFF真彩色合成def write_rgb_geotiff(filepath, r_array, g_array, b_array, geotransform, projection): r/g/b_array: 均为 2D numpy 数组形状相同 driver gdal.GetDriverByName(GTiff) ds driver.Create(filepath, r_array.shape[1], r_array.shape[0], 3, gdal.GDT_UInt16) ds.SetGeoTransform(geotransform) ds.SetProjection(projection) # 写入 R、G、B 波段1-indexed ds.GetRasterBand(1).WriteArray(r_array) ds.GetRasterBand(2).WriteArray(g_array) ds.GetRasterBand(3).WriteArray(b_array) # 设置波段描述可选但利于后续识别 ds.GetRasterBand(1).SetDescription(Red Band) ds.GetRasterBand(2).SetDescription(Green Band) ds.GetRasterBand(3).SetDescription(Blue Band) ds.FlushCache() del ds # 合成 Sentinel-2 真彩色B04Red, B03Green, B02Blue write_rgb_geotiff(sentinel2_truecolor.tif, b04, b03, b02, gt, proj)参数说明Create()第四参数是数据类型gdal.GDT_UInt16适配 12bit 遥感数据若用Float32需确保array是 float 类型否则静默截断。4. 常见故障排查为什么读出来全是 0为什么写入后 QGIS 打不开生产环境中80% 的 GDAL 问题源于元数据缺失、类型不匹配或路径权限。4.1 读取全零或异常值的 5 个检查点检查项命令/代码期望结果问题定位1. 文件是否被其他进程锁定lsof | grep your_file.tif(Linux/macOS) 或资源监视器 (Windows)无输出文件被 ENVI/QGIS 占用GDAL 只读模式失败2. 波段数据类型是否匹配print(band.DataType)1(Byte),2(UInt16),6(Float32)若误用ReadAsArray()读UInt16到int8数组高位被截断为 03. NoData 值是否被忽略print(band.GetNoDataValue())如0,65535,-9999读取后需array[array nodata] np.nan4. 图像压缩是否不支持print(ds.GetMetadata(IMAGE_STRUCTURE))COMPRESSION: LZW或DEFLATEGDAL 3.4 支持 LZW/DEFLATE旧版需--config GDAL_PAM_ENABLED NO5. 内存映射是否启用gdal.Info(file.tif, formatjson)[metadata][IMAGE_STRUCTURE]INTERLEAVE: BAND或PIXELBAND交错存储更省内存PIXEL适合 OpenCV 处理4.1.1 快速诊断脚本保存为gdal_check.pyfrom osgeo import gdal import sys def check_tiff(filepath): ds gdal.Open(filepath) if ds is None: print(f❌ 打开失败{filepath}) return print(f✅ 文件打开成功) print(f 驱动: {ds.GetDriver().ShortName}) print(f 尺寸: {ds.RasterXSize}x{ds.RasterYSize}x{ds.RasterCount}) gt ds.GetGeoTransform() print(f GeoTransform: {gt[:4]}... (前4位)) proj ds.GetProjection() print(f 投影长度: {len(proj)} 字符) for i in range(1, min(3, ds.RasterCount 1)): b ds.GetRasterBand(i) dt gdal.GetDataTypeName(b.DataType) ndv b.GetNoDataValue() print(f 波段{i}: {dt}, NoData{ndv}) # 读取左上角 3x3 像元检查值 win b.ReadAsArray(0, 0, 3, 3) print(f 左上角值: {win.flatten()}) if __name__ __main__: check_tiff(sys.argv[1])运行python gdal_check.py landsat8_b5.tif输出清晰显示每层状态比gdalinfo更聚焦遥感处理痛点。4.2 写入后 QGIS/GIS 软件打不开的硬核修复QGIS 报错 “Invalid layer” 或 “Projection not recognized” 通常因 WKT 投影字符串不标准。4.2.1 强制标准化投影解决 EPSG 代号丢失from osgeo import osr def standardize_projection(wkt_str): 将任意 WKT 投影转为标准 EPSG 代号格式 srs osr.SpatialReference() srs.ImportFromWkt(wkt_str) # 尝试导出为 EPSG 代号如存在 epsg_code srs.AutoIdentifyEPSG() if epsg_code 0: return wkt_str # 无法识别返回原字符串 else: srs_out osr.SpatialReference() srs_out.ImportFromEPSG(int(srs.GetAuthorityCode(None))) return srs_out.ExportToWkt() # 修复后写入 standard_proj standardize_projection(proj) write_geotiff(fixed.tif, array, gt, standard_proj)4.2.2 添加地理参考当 GeoTransform 为 (0,1,0,0,0,1) 时def fix_identity_geotransform(ds, lon_min, lat_min, lon_max, lat_max): 用地理范围反推仿射变换适用于无地理参考的 TIFF xsize, ysize ds.RasterXSize, ds.RasterYSize xres (lon_max - lon_min) / xsize yres (lat_max - lat_min) / ysize # 左上角坐标 gt (lon_min, xres, 0, lat_max, 0, -yres) ds.SetGeoTransform(gt) # 设置 WGS84 投影 srs osr.SpatialReference() srs.ImportFromEPSG(4326) ds.SetProjection(srs.ExportToWkt()) return ds # 示例给一张纯 RGB 图添加 WGS84 地理参考 ds gdal.Open(rgb.jpg, gdal.GA_Update) ds fix_identity_geotransform(ds, 116.2, 39.7, 116.5, 40.0) ds.FlushCache() del ds注意gdal.GA_Update模式打开才能写入元数据。若原文件只读需先shutil.copy()到新路径再操作。5. 进阶技巧用 GDAL 处理国产高分影像的坐标偏移与波段顺序国产高分系列GF-1/2/6/7影像常含两个特殊问题1WGS84 坐标系下存在百米级偏移2波段顺序与国际标准不一致如 GF-1 PMS 的蓝绿红近红外顺序是 B2,B3,B4,B5而非 B1,B2,B3,B4。GDAL 提供原生接口修正。5.1 用 GCPs地面控制点校正坐标偏移高分影像的.xml附带 GCP 列表GDAL 可直接读取并生成 RPC 模型。5.1.1 从 XML 提取 GCPs 并写入 GeoTIFFimport xml.etree.ElementTree as ET def parse_gf_gcps(xml_path): 解析高分 XML 中的 GCPs tree ET.parse(xml_path) root tree.getroot() gcps [] for gcp in root.findall(.//GCP): line float(gcp.find(Line).text) samp float(gcp.find(Sample).text) lat float(gcp.find(Lat).text) lon float(gcp.find(Lon).text) height float(gcp.find(Height).text) if gcp.find(Height) is not None else 0 gcps.append(gdal.GCP(lon, lat, height, samp, line)) return gcps # 将 GCPs 写入 TIFF替代 GeoTransform ds gdal.Open(gf1_pms1.tif, gdal.GA_Update) gcps parse_gf_gcps(GF1_PMS1_E116.3_N39.8_20230501_L1A0000111111.xml) ds.SetGCPs(gcps, WGS84) # 生成 .aux.xml 辅助文件QGIS 识别 GCP ds.BuildOverviews(NEAREST, [2,4,8]) ds.FlushCache() del ds效果QGIS 加载后自动启用 GCP 校正消除系统性偏移。BuildOverviews生成金字塔加速大图浏览。5.2 自动识别高分波段并重排顺序高分数据常以B*命名波段但 GDAL 默认按文件内顺序读取。需解析元数据确定物理波长。5.2.1 读取波段中心波长并排序def get_band_wavelengths(ds): 从 GDAL 元数据或 XML 解析各波段中心波长nm # 方法1检查 GDAL 元数据 meta ds.GetMetadata() if BAND_WAVELENGTH in meta: return [float(x) for x in meta[BAND_WAVELENGTH].split(,)] # 方法2解析附带 XML高分标准 import glob xml_files glob.glob(ds.GetDescription().replace(.tif, *.xml)) if xml_files: tree ET.parse(xml_files[0]) wavelengths [] for band in tree.findall(.//Band): wl band.find(CenterWavelength) if wl is not None: wavelengths.append(float(wl.text)) return wavelengths # 默认波长GF-1 PMS return [450, 520, 630, 830] # B2,B3,B4,B5 def reorder_bands_to_rgb(ds, target_order[2, 1, 0]): # B4(R), B3(G), B2(B) 按波长重排波段生成真彩色数组 wavelengths get_band_wavelengths(ds) # 按波长升序排序索引蓝绿红近红外 sorted_idx sorted(range(len(wavelengths)), keylambda i: wavelengths[i]) # 映射到目标顺序取第2个红、第1个绿、第0个蓝 rgb_idx [sorted_idx[i] for i in target_order] data np.stack([ ds.GetRasterBand(i1).ReadAsArray().astype(np.float32) for i in rgb_idx ], axis0) return data # 使用 ds gdal.Open(GF1_PMS1.tif) rgb_data reorder_bands_to_rgb(ds) # 自动匹配 B4,B3,B2 write_rgb_geotiff(gf1_truecolor.tif, *rgb_data, ds.GetGeoTransform(), ds.GetProjection())逻辑说明target_order[2,1,0]表示取波长第三长红、第二长绿、第一长蓝的波段。sorted_idx是按波长从小到大排列的原始波段索引rgb_data保证输出符合人眼感知。至此你已掌握用 Python GDAL 处理 GeoTIFF 遥感影像的全链路从避坑安装、分块读取、地理裁剪、带参写入到国产数据偏移校正与波段重排。下一步可延伸至gdal.Warp投影转换、gdal.Translate格式转换、或接入rasterio做更 Pythonic 的操作——但核心原理始终是GeoTransform、Projection、RasterBand这三根支柱。本文还有配套的精品资源点击获取
