前言

最近使用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

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

python使用gdal读取tif经纬度
python使用gdal读取tif经纬度前言一.tif是什么? 二.使用gdal读取经纬度 1.引入库 2.读取坐标总结前言博主作为一个GIS开发者,开发过程中不免遇到一些处理tif的问题和 ...
python安装gdal包_python安装gdal的两种方法
1.不用手动下载文件,直接执行以下命令即可 conda install gdal 2.首先,下载gdal的whl文件链接, 官网下载比较慢,GDAL-2.2.4-cp27-cp27m-win_am ...
Python扩展库scipy.misc中图像转换成pillow图像
众所周知,在数字图像处理领域中有很多基准测试图像,这些图像用来作为科研人员PK自己的算法时的参考,给大家提供一个公平的样本,针对同一个问题进行处理时,可以用这些基准图像做实验,比较常见的应该就是len ...
python中的scipy基础知识_python中SciPy是什么？
python中Numpy常用于计算二维数组计算,而python的另一个库SciPy库与Numpy有着密切的关系,是需要通过Numpy为基础,同时也是通过Numpy数据来操控科学计算.常见的是插值运算. ...
python使用GDAL/OGR/OSR时设置GDAL_DATA环境变量
版权声明:转载请注明作者(独孤尚良dugushangliang)出处: https://blog.csdn.net/dugushangliang/article/details/89377625 右键 ...
Python库之Scipy库的简介、安装详细
目录 Scipy库的简介 Scipy库的安装 Scipy库的简介 Scipy高级科学计算库:和Numpy联系很密切,Scipy一般都是操控Numpy数组来进行科学计算.统计分析,所以可以说是基于Num ...
Python使用GDAL矢量裁剪栅格，设置背景值为空白（已解决）
一.使用gdal.Warp gdalwarp 实用程序是一种图像拼接.重投影和扭曲实用程序.该程序可以重新投影到任何支持的投影.如果图像是带有控制信息的"原始"图像,也可以存储原始 ...
python中使用scipy.integrate求积分、二重积分、三重积分
python中使用scipy.integrate求积分.二重积分.三重积分代码如下: import numpy as np from scipy.integrate import quad, tpl ...
python统计函数库scipy.stats的用法1/3
背景总结统计工作中几个常用用法在python统计函数库scipy.stats的使用范例. 正态分布以正态分布的常见需求为例了解scipy.stats的基本使用方法. 生成服从指定分布的随机数 no ...

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