根据经纬度获取海拔高度

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