Python使用GDAL读取GeoTiff时查找指定经纬度对应最近像素的方法

核心思路 #

• 首先将输入的WGS84经纬度坐标，转换为GeoTiff本身使用的投影坐标系坐标
• 通过GDAL内置的逆地理变换能力，将投影坐标转换为栅格对应的像素行列号
• 判断行列号是否在栅格有效范围内，返回对应标识与像素值

完整实现代码 #

from osgeo import gdal, osr
import numpy as np

def get_pixel_from_lonlat(ds, lon, lat):
    # 读取栅格基本信息
    width = ds.RasterXSize
    height = ds.RasterYSize
    gt = ds.GetGeoTransform()
    proj = ds.GetProjection()
    
    # 1. 坐标转换：WGS84经纬度转栅格投影坐标
    # 定义源坐标系（WGS84 EPSG:4326）
    src_srs = osr.SpatialReference()
    src_srs.ImportFromEPSG(4326)
    # 定义目标坐标系（栅格自身投影）
    dst_srs = osr.SpatialReference()
    dst_srs.ImportFromWkt(proj)
    # 适配GDAL3+的轴顺序规则，强制按传统GIS的经度在前、纬度在后顺序处理
    if int(gdal.VersionInfo()[0]) >= 3:
        src_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
        dst_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    transform = osr.CoordinateTransformation(src_srs, dst_srs)
    # 执行转换，返回格式为(投影x, 投影y, 高程)
    x, y, _ = transform.TransformPoint(lon, lat)
    
    # 2. 逆地理变换，投影坐标转像素行列号
    success, inv_gt = gdal.InvGeoTransform(gt)
    if not success:
        return False, None, None, None
    px = int(round(inv_gt[0] + inv_gt[1] * x + inv_gt[2] * y))
    py = int(round(inv_gt[3] + inv_gt[4] * x + inv_gt[5] * y))
    
    # 3. 范围校验
    if px < 0 or px >= width or py < 0 or py >= height:
        return False, px, py, None
    
    # 4. 仅读取目标像素值，避免加载全量数据，大文件场景性能更优
    data = ds.ReadAsArray(px, py, 1, 1)
    return True, px, py, data.squeeze()

# 使用示例
if __name__ == "__main__":
    ds = gdal.Open('foo.tiff')
    # 替换为你要查询的经纬度
    target_lon = -95
    target_lat = 40
    is_valid, px, py, pixel_value = get_pixel_from_lonlat(ds, target_lon, target_lat)
    if is_valid:
        print(f"坐标在栅格范围内，像素行列号：({px}P, {py}L)，像素值：{pixel_value}")
    else:
        print("目标坐标超出栅格覆盖范围")

相关说明 #

• 该实现逻辑和gdallocationinfo -wgs84命令完全一致，依赖GDAL内置能力完成坐标转换与像素定位，不需要额外生成经纬度网格或暴力搜索，性能远优于后备方案
• 兼容带旋转参数的GeoTiff，适配绝大多数栅格文件场景
• 若你的GDAL为旧版本，导入时报错可将导入语句改为：

  import gdal
  import osr
  

内容的提问来源于stack exchange，提问作者Grant Petty