利用经纬度查询海拔高度

1、准备 安装python包

apt-get update

apt-get install -y python3-gdal

pip install gdal

pip install pyproj

2、下载高程tif数据

全国各省12.5m高程数据-免费下载-资源下载-数字地球开放平台

注册下，下载他的网盘下载工具，然后选取要使用的区域下载。我下载的宁夏的。

3、python代码 

以下的代码是专门针对以上的高程数据写的代码。不同的高程数据可能需要不同的投影计算方法。

使用Pyproj库进行投影转换:

        通过pyproj.Proj方法，确认源投影（EPSG:4326）和目标投影（EPSG:32648）。

        将输入的经纬度数据从WGS 84转换为UTM 48N坐标系统。

import gdal
import osr
from pyproj import Proj, transform

def print_geo_transform_info(transform):
    print(f"Geotransform:")
    print(f"  Origin: ({transform[0]}, {transform[3]})")
    print(f"  Pixel Size: ({transform[1]}, {transform[5]})")

def get_elevation(tif_path, lon, lat):
    # 打开 GeoTIFF 文件
    dataset = gdal.Open(tif_path)
    if dataset is None:
        raise FileNotFoundError(f"Cannot open file: {tif_path}")

    # 打印GeoTIFF文件信息
    print("GeoTIFF File Information:")
    print(f"Raster size: {dataset.RasterXSize} x {dataset.RasterYSize}")
    print(f"Number of bands: {dataset.RasterCount}")

    # 获取地理变换信息
    transform_info = dataset.GetGeoTransform()
    if transform_info is None:
        raise ValueError("Failed to get geotransform from the dataset.")
    print_geo_transform_info(transform_info)

    # 获取投影信息
    proj_info = dataset.GetProjection()
    if proj_info is None:
        raise ValueError("Failed to get projection from the dataset.")
    print(f"Projection: {proj_info}")

    # 获取栅格波段
    band = dataset.GetRasterBand(1)
    if band is None:
        raise ValueError("Failed to get raster band from the dataset.")
    
    # 创建投影转换对象
    in_proj = Proj(init='epsg:4326')  # WGS 84
    out_proj = Proj(init='epsg:32648')  # UTM zone 48N

    try:
        # 将经纬度坐标转换为UTM坐标
        srcX, srcY = transform(in_proj, out_proj, lon, lat)
        print(f"Transformed Coordinates: X={srcX}, Y={srcY}")

        # 检查是否在GeoTIFF数据范围内
        x_min = transform_info[0]
        x_max = x_min + transform_info[1] * dataset.RasterXSize
        y_max = transform_info[3]
        y_min = y_max + transform_info[5] * dataset.RasterYSize

        print(f"GeoTIFF X range: {x_min} to {x_max}")
        print(f"GeoTIFF Y range: {y_min} to {y_max}")

        if not (x_min <= srcX <= x_max and y_min <= srcY <= y_max):
            raise ValueError("The given coordinates are outside the GeoTIFF boundaries.")

        # 将栅格坐标转换为像素坐标
        px = int((srcX - transform_info[0]) / transform_info[1])
        py = int((srcY - transform_info[3]) / transform_info[5])

        # 获取像素值（高程值）
        elevation = band.ReadAsArray(px, py, 1, 1)[0][0]
        return elevation
    except Exception as e:
        print(f"Coordinate transformation failed: {e}")
        return None

if __name__ == "__main__":
    # GeoTIFF 文件路径
    tif_path = "/usr/local/yxq/map/ningxia.tif"
    # 经纬度坐标
    lon = 106.0836457924024
    lat = 36.218175877348614

    try:
        # 获取高程数据
        elevation = get_elevation(tif_path, lon, lat)
        if elevation is not None:
            print(f"经度: {lon}, 纬度: {lat}, 海拔高度: {elevation} 米")
        else:
            print(f"Cannot determine elevation for coordinates 经度: {lon}, 纬度: {lat}")
    except Exception as e:
        print(f"Error: {e}")

作者：yixiuquan