（Python）使用Gdal+Scipy获得Dem的经纬度的高程值（双线性和三次样条内插）

前言

最近使用python进行一些的遥感影像处理，发现摄影测量方面相比于CV领域资料非常匮乏，特此记录一些小工具的实现，也为遥感领域的开源做出自己的一些贡献，水平有限，能用而已。

基本原理

Dem内插是指根据分布在内插点周围的已知参考点的高程值求出未知点的高程值，由于所采集的原始数据排列一般是不规则的，为了获取规则的 DEM，内插是必不可少的重要步骤。 它是 DEM的核心问题，贯穿于 DEM 的生产、质量控制、精度评定、分析应用的各个环节。

双线性内插

三次样条内插

代码实现

from osgeo import gdal
import numpy as np
import math
from scipy import interpolate
#func：lon，lat到像素坐标
def geo2imagexy(trans, x, y):
    a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
    b = np.array([x - trans[0], y - trans[3]])
    return np.linalg.solve(a, b)  # 使用numpy的linalg.solve进行二元一次方程的求解

#读入dem
dem = gdal.Open('D:\\DEM_EGM96.tif')
dem_geotrans = dem.GetGeoTransform()  # 地理信息仿射矩阵
band = dem.GetRasterBand(1)
elevation = band.ReadAsArray()#灰度值读进来
test_geo = [112.961403969,34.452284651]
y,x= geo2imagexy(dem_geotrans, test_geo [0], test_geo [1])
test_point = [x,y]
# #双线性插值内插出高程
# #取最近的四个点
# x_begin, x_end = math.floor(test_point[0]), math.ceil(test_point[0])
# y_begin, y_end = math.floor(test_point[1]), math.ceil(test_point[1])
# z = elevation[x_begin:x_end + 1, y_begin:y_end + 1]
# func = interpolate.interp2d([x_begin, x_end],[y_begin, y_end], z, kind='linear')
# value = func(test_point[0], test_point[1])

# 三次样条二维插值内插出高程
#取周围 4*4的点
x_begin, x_end = math.floor(test_point[0])-1, math.ceil(test_point[0])+1
y_begin, y_end = math.floor(test_point[1])-1, math.ceil(test_point[1])+1
z = elevation[x_begin:x_end + 1, y_begin:y_end + 1]
xnew=np.arange(x_begin, x_end + 1)
ynew=np.arange(y_begin, y_end + 1)
func = interpolate.interp2d(xnew,ynew, z, kind='cubic')
value = func(test_point[0], test_point[1])
print(value)

后记

代码实现比较简陋，也没有做边界判断，就是提供一个实现思路。Envi上测试了几个点，基本正确，如果问题，一起讨论，后续还有分块匹配，RPC相关的读取、正投反投、三角化的相关实现，敬请期待。
如有交流，个人联系方式：WHUwsd1995