Python GDAL 根据tif文件 计算高程数据，两点间的水平距离

知道 高程数据，两点间的水平距离就可以计算坡度了

from osgeo import gdal
from osgeo import osr
import numpy as np
import sys
import math
import datetime
from flask import Flask, request
from pyproj import Transformer



# 读取tif文件
def readHuNanTif():

    gdal.AllRegister()
    dataset = gdal.Open("D:/hunan.tif")

    rasterCount = dataset.RasterCount
    # print(dataset.RasterCount)  # 波段数

    cols = dataset.RasterXSize  # 图像长度
    rows = (dataset.RasterYSize)  # 图像宽度

    xoffset = cols
    yoffset = rows
    print('图像宽', xoffset)
    print('图像高', yoffset)


    # 读取光栅数据
    band = dataset.GetRasterBand(1)

    # 高程数据
    print('开始读取', datetime.datetime.now())
    ele = band.ReadAsArray()
    print('读取结束', datetime.datetime.now())
    return ele


# 根据经纬度读取高程信息
def readHeight(lon, lat):

    nrows, ncols = elevation.shape

    x0, dx, dxdy, y0, dydx, dy = dataset.GetGeoTransform()

    # 根据经纬度计算行列号，dx=dy为分辨率，不相等的时候(y0-latitude)/dx改为(y0-latitude)/-dy
    new_ncols, new_nrows = int((y0 - lat) / dx), int((lon - x0) / dx)

    # 根据行列号读取并打印输出指定坐标点高程
    print('读高程数据', datetime.datetime.now())
    height = elevation[new_ncols][new_nrows]
    print('高程结束中', datetime.datetime.now())
    print('高程数据', height)

    return height


# 计算两点的水平距离
def calculateHorizontalDistance(hb_lon, hb_lat, ha_lon, ha_lat):
    point_trans_b = axis_conversion(hb_lon, hb_lat)
    print('投影坐标', point_trans_b)

    point_trans_a = axis_conversion(ha_lon, ha_lat)
    print('投影坐标', point_trans_a)

    # 根号    (xa - xb) 平方 + (ya - yb) 平方
    xa = point_trans_a[0]
    xb = point_trans_b[0]
    ya = point_trans_a[1]
    yb = point_trans_b[1]

    dis = math.sqrt((xa - xb) * (xa - xb) + (ya - yb) * (ya - yb))

    return dis


# 经纬度转 投影坐标
def axis_conversion(lon, lat):
    transformer = Transformer.from_crs("epsg:4326", "epsg:3857")
    x3, y3 = transformer.transform(lat, lon)
    return (x3, y3)



readHuNanTif()