Python+GDAL生成数字表面模型全流程（持续更新中）

   在卫星摄影测量中，在没有控制点的情况下生成DEM，可以利用卫星影像附带的RPC文件来辅助生成DEM。其主要分为

1、核线影像生成（各种各样的方法）

2、密集匹配（一般采用SGM）

3、RPC模型平差

在第一步的核线影像生成中，在不知道卫星参数的情况下，一般利用RPC辅助生成核线影像。

首先要读取卫星原始数据中的RPC参数，本次利用python+GDAL读取了卫星影像的RPC参数。

具体代码如下：

import os
import re
import numpy as np
import string
from osgeo import gdal
import time

# the class of RPC, including the par and function of rpc
class RPC():
    lineNumCoef = []
    lineDenCoef = []
    sampNumCoef = []
    sampDenCoef = []
    errBias = 0
    errRand = 0
    lineOffset = 0
    sampOffset = 0
    latOffset = 0
    longOffset = 0
    heightOffset = 0
    lineScale = 0
    sampScale = 0
    latScale = 0
    longScale = 0
    heightScale = 0

def pre_process_RPC(rpc):
    sample = RPC()
    sample.errBias = rpc['ERR_BIAS']
    sample.errRand = rpc['ERR_RAND']
    sample.lineOffset = rpc['LINE_OFF']
    sample.sampOffset = rpc['SAMP_OFF']
    sample.latOffset = rpc['LAT_OFF']
    sample.longOffset = rpc['LONG_OFF']
    sample.heightOffset = rpc['HEIGHT_OFF']
    sample.lineScale = rpc['LINE_SCALE']
    sample.sampScale = rpc['SAMP_SCALE']
    sample.latScale = rpc['LAT_SCALE']
    sample.longScale = rpc['LONG_SCALE']
    sample.heightScale = rpc['HEIGHT_SCALE']

    string1 = rpc['SAMP_DEN_COEFF'].split()
    string2 = rpc['SAMP_NUM_COEFF'].split()
    string3 = rpc['LINE_DEN_COEFF'].split()
    string4 = rpc['LINE_NUM_COEFF'].split()

    for i in range(20):
        sample.sampDenCoef.append(float(string1[i]))
        sample.sampNumCoef.append(float(string2[i]))
        sample.lineDenCoef.append(float(string3[i]))
        sample.lineNumCoef.append(float(string4[i]))

    return sample

input = "**************"    # path of the rpc file
dataset = gdal.Open(input, gdal.GA_Update)
rpc = dataset.GetMetadata("RPC")
result = pre_process_RPC(rpc)

在利用投影轨迹法计算核线轨迹的时候需要用到RPC模型，具体公式如下

RPC 模型正解公式

 具体的Python代码如下：

# Get l and s by RPC model
def RPC_Positive(lat, long, height, rpc):
    U = (lat - rpc.latOffset) / rpc.latScale
    V = (lon - rpc.lonOffset) / rpc.longScale
    W = (height - rpc.heightOffset) / rpc.heightScale

    u = np.mat('1, V, U, W, V*U, V*W, U*W, V*V, U*U, W*W,\
     U*V*W, V*V*V, V*U*U, V*W*W, V*V*U,\
     U*U*U, U*W*W, V*V*W, U*U*W, W*W*W')
    a = np.mat(rpc.lineNumCoef)
    b = np.mat(rpc.lineDenCoef)
    c = np.mat(rpc.sampleNumCoef)
    d = np.mat(rpc.sampleDenCoef)

    X = np.matmul(a, u.T) / np.matmul(b, U.t)
    Y = np.matmul(c, u.T) / np.matmul(d, u.T)

    l = rpc.lineOffset + rpc.lineScale * X
    s = rpc.sampleOffset + rpc.sampleScale * Y

    return l, s

下面还需要RPC反解算公式（由行列号计算经纬度）

直接放代码：

def RPC_Negative_inverse(samp, line , height, rpc):

    Xf = np.vstack([line, samp]).T

    long = -line ** 0
    lat = -samp ** 0

    EPS = 2
    x0 = cal_rfm(rpc.lineNumCoef, rpc.lineDenCoef, lat, long, height)
    y0 = cal_rfm(rpc.sampNumCoef, rpc.sampDenCoef, lat, long, height)
    x1 = cal_rfm(rpc.lineNumCoef, rpc.lineDenCoef, lat, long + EPS, height)
    y1 = cal_rfm(rpc.sampNumCoef, rpc.sampDenCoef, lat, long + EPS, height)
    x2 = cal_rfm(rpc.lineNumCoef, rpc.lineDenCoef, lat + EPS, long, height)
    y2 = cal_rfm(rpc.sampNumCoef, rpc.sampDenCoef, lat + EPS, long, height)

    n = 0
    while not np.all((x0 - line) ** 2 + (y0 - samp) ** 2 < 1e-20):

        if n > 1000:
            print("error,exceed 1000 times RPC_Negative")
            pass
        X0 = np.vstack([x0, y0]).T
        X1 = np.vstack([x1, y1]).T
        X2 = np.vstack([x2, y2]).T
        e1 = X1 - X0
        e2 = X2 - X0
        u = Xf - X0

        num = np.sum(np.multiply(u, e1), axis=1)
        den = np.sum(np.multiply(e1, e1), axis=1)
        a1 = np.divide(num, den).squeeze()
        num = np.sum(np.multiply(u, e2), axis=1)
        den = np.sum(np.multiply(e2, e2), axis=1)
        a2 = np.divide(num, den).squeeze()

        long = long + a1 * EPS
        lat = lat + a2 * EPS

        EPS = .1
        x0 = cal_rfm(rpc.lineNumCoef, rpc.lineDenCoef, lat, long, height)
        y0 = cal_rfm(rpc.sampNumCoef, rpc.sampDenCoef, lat, long, height)
        x1 = cal_rfm(rpc.lineNumCoef, rpc.lineDenCoef, lat, long + EPS, height)
        y1 = cal_rfm(rpc.sampNumCoef, rpc.sampDenCoef, lat, long + EPS, height)
        x2 = cal_rfm(rpc.lineNumCoef, rpc.lineDenCoef, lat + EPS, long, height)
        y2 = cal_rfm(rpc.sampNumCoef, rpc.sampDenCoef, lat + EPS, long, height)

        n += 1
    return long, lat


def RPC_Negative(line, samp, height, rpc):
    line = (line - rpc.lineOffset) / rpc.lineScale
    samp = (samp - rpc.sampOffset) / rpc.sampScale
    height = (height - rpc.heightOffset) / rpc.heightScale

    long, lat = RPC_Negative_inverse(samp, line, height, rpc)

    long = long * rpc.longScale + rpc.longOffset
    lat = lat * rpc.latScale + rpc.latOffset

    return float(long), float(lat)

最近发现一个事情，GDAL库里面提供了rpc正反算，我自己编写似乎没太大意义。