（Python）卫星RPC有理多项式模型读取与正反投影坐标计算原理与实现（Python）卫星RPC有理多项式模型读取与正反投影坐标计算原理与实现
⽂章⽬录
摘要
在航天摄影测量领域，RPC模型是⼀种通⽤的卫星⼏何定位模型。作为基础构件，RPC模型是后续前⽅交会（三⾓化）、区域⽹平差、DSM⽣成的基础。但是，笔者在实际使⽤中发现⽹上的关于RPC模型开源代码良莠不齐，有的介绍详细但没有代码，难以找到⼀组⽐较易⽤的RPC模型库。经过寻找测试，笔者发现s2p作者开源的rpm具有健全的功能与清晰的代码，笔者将其稍微修改与整理，并添加了注释，使其变为依赖少、可应⽤⼤部分国产卫星RPC处理的单⽂件功能库，功能主要包括：1.读取RPC模型；2.RPC模型正投影，即地理坐标计算像素坐标；3.RPC模型逆投影，即像素坐标与⾼程计算经纬度。后续也将基于这个简单的⼩库，做⼀些应⽤的演⽰。
本⽂⾸先介绍了RPC⼏何定位模型的基本知识，然后提供了RPC模型代码的实现并进⾏了简单的使⽤⽰范，最后评估的该代码的精度与性能。
RPC⼏何定位模型介绍祷告的意思
卫星遥感影像在成像过程中由于受到诸多复杂因素的影响，使各像点产⽣了不同程度的⼏何变形。建⽴遥感影像⼏何定位模型可以正确地描述每⼀个像点坐标与其对应地⾯点物⽅坐标间的严格⼏何关系，以便对原始影像进⾏⾼精度的⼏何纠正及对地⽬标定位，从⽽实现由⼆维影像反演实地表⾯的平⾯或空间位置，以满⾜各种遥感应⽤的需求。
现今已建⽴了各种传感器模型，来描述地⾯空间坐标与相应像点坐标的数学关系.在传统的摄影测量领域，应⽤较多的是物理模型.这种模型的处理技术已趋向成熟，定位精度⽐较⾼.但由于物理传感器模型涉及传感器物理结构，成像⽅式及各种成像参数.为了保护技术秘密,⼀些公司只使⽤有理函数模型(RFM:Rational Function Model)，并向⽤户提供RFM的参数――有理多项式系数(RPC:Rational Polynomial Coefficients).
RPC可以将其理解为⼀个带畸变参数的相机模型，其描述的是从三维的地理坐标到⼆维的卫星影像坐标之间的转换关系，⼀般称之为从物⽅到像⽅，我们可以理解为3D到2D，这也被称之为正投影，具体为已知lon，lat，h，求得像素坐标s，l。⽽反向投影的则是从影像坐标系到地理坐标系，即从像⽅到物⽅，这⾥要注意的是，这仍然是⼀个3D到2D的关系，因为2D是⽆法升维到3D的，具体为已经s，l，h，求得lon与lat。很多同学会在这个地⽅犯迷糊，注意已知什么，求什么，就可以清楚把握。
RPC模型正投影表达式如下(截图来⾃：杨博，王密，⽪英冬．仅⽤虚拟控制点的超⼤区域⽆控制区域⽹平差［Ｊ］．测绘学报）
[外链图⽚转存失败，源站可能有防盗链机制，建议将图⽚保存下来直接上传(img-nYRWrwR1-1637063158167)(C:\Urs\卫少东
\Documents\BaiduNetdiskWorkspace\个⼈⽂档\博客写作\国产卫星RPC有理多项式模型读取与正反投影坐标计算
\body.asts\image-20211116174102785.png)]
逆投影（反投影）则是将地理坐标当做待求解参数，进⾏平差迭代求解，具体原理各⼤摄影教科书都有，此处不再赘述。
RPC模型库代码实现
Rpc模型库
#命名为：geo_rpc.py
"""
RPC model parrs, localization, and projection
"""
import numpy as np
from osgeo import gdal
#最⼤迭代次数超过则报错
class MaxLocalizationIterationsError(Exception):
"""
Custom rpcm Exception.
"""
pass
def apply_poly(poly, x, y, z):
"""
Evaluates a 3-variables polynom of degree 3 on a triplet of numbers.
将三次多项式的统⼀模式构建为⼀个单独的函数
Args:
表达方式的作用poly: list of the 20 coefficients of the 3-variate degree 3 polynom,
ordered following the RPC convention.
x, y, z: triplet of floats. They may be numpy arrays of same length.
Returns:
the value(s) of the polynom on the input point(s).
"""
out =0
out += poly[0]
out += poly[1]*y + poly[2]*x + poly[3]*z
out += poly[4]*y*x + poly[5]*y*z +poly[6]*x*z
out += poly[7]*y*y + poly[8]*x*x + poly[9]*z*z
out += poly[10]*x*y*z
out += poly[11]*y*y*y
out += poly[12]*y*x*x + poly[13]*y*z*z + poly[14]*y*y*x
out += poly[15]*x*x*x
out += poly[16]*x*z*z + poly[17]*y*y*z + poly[18]*x*x*z
out += poly[19]*z*z*z
return out
def apply_rfm(num, den, x, y, z):
"""
Evaluates a Rational Function Model (rfm), on a triplet of numbers.
执⾏20个参数的分⼦和20个参数的除法
Args:
num: list of the 20 coefficients of the numerator
den: list of the 20 coefficients of the denominator
All the coefficients are ordered following the RPC convention.
x, y, z: triplet of floats. They may be numpy arrays of same length.
Returns:
the value(s) of the rfm on the input point(s).
"""
return apply_poly(num, x, y, z)/ apply_poly(den, x, y, z)
def rpc_from_geotiff(geotiff_path):
"""
Read the RPC coefficients from a GeoTIFF file and return an RPCModel object.该函数返回影像的Gdal格式的影像和RPCmodel
Args:
geotiff_path (str): path or url to a GeoTIFF file
Returns:
instance of the rpc_model.RPCModel class
"""
coagulation
# with rasterio.open(geotiff_path, 'r') as src:
#
datat = gdal.Open(geotiff_path, gdal.GA_ReadOnly)
rpc_dict = datat.GetMetadata("RPC")
# 同时返回影像与rpc
return datat, RPCModel(rpc_dict,'geotiff')
def read_rpc_file(rpc_file):
"""
Read RPC from a RPC_txt file and return a RPCmodel
从TXT中直接单独读取RPC模型
Args:
rpc_file: RPC sidecar file path
Returns:
chopsticks
dictionary read from the RPC file, or an empty dict if fail
"""
with open(rpc_file)as f:
with open(rpc_file)as f:
rpc_content = f.read()
rpc = read_rpc(rpc_content)
return RPCModel(rpc)
def read_rpc(rpc_content):
"""
Read RPC file assuming the ikonos format
解析RPC参数
Args:
rpc_content: content of RPC sidecar file path read as a string
Returns:
dictionary read from the RPC file
"""
import re
lines = rpc_content.split('\n')
ceo是什么dict_rpc ={}
for l in lines:
ll = l.split()
if len(ll)>1:
k = re.sub(r"[^a-zA-Z0-9_]","",ll[0])
dict_rpc[k]= ll[1]
def par_coeff(dic, prefix, indices):
""" helper function"""
return' '.join([dic["%s_%s"%(prefix,str(x))]for x in indices])
dict_rpc['SAMP_NUM_COEFF']= par_coeff(dict_rpc,"SAMP_NUM_COEFF",range(1,21)) dict_rpc['SAMP_DEN_COEFF']= par_coeff(dict_rpc,"SAMP_DEN_COEFF",range(1,21)) dict_rpc['LINE_NUM_COEFF']= par_coeff(dict_rpc,"LINE_NUM_COEFF",range(1,21)) dict_rpc['LINE_DEN_COEFF']= par_coeff(dict_rpc,"LINE_DEN_COEFF",range(1,21))
return dict_rpc
class RPCModel:
def__init__(lf, d, dict_format="geotiff"):
"""
Args:
d (dict): dictionary read from a geotiff fil
e with
rasterio.open('/path/to/file.tiff', 'r').tags(ns='RPC'),
or from the .__dict__ of an RPCModel object.
dict_format (str): format of the dictionary pasd in `d`.
Either "geotiff" if read from the tags of a geotiff file,
or "rpcm" if read from the .__dict__ of an RPCModel object.
"""
if dict_format =="geotiff":韵文
lf.lat_offt =float(d['LAT_OFF'][0:d['LAT_OFF'].rfind(' ')])
lf.lon_offt =float(d['LONG_OFF'][0:d['LONG_OFF'].rfind(' ')])
lf.alt_offt =float(d['HEIGHT_OFF'][0:d['HEIGHT_OFF'].rfind(' ')])
lf.lat_scale =float(d['LAT_SCALE'][0:d['LAT_SCALE'].rfind(' ')])
lf.lon_scale =float(d['LONG_SCALE'][0:d['LONG_SCALE'].rfind(' ')])
lf.alt_scale =float(d['HEIGHT_SCALE'][0:d['HEIGHT_SCALE'].rfind(' ')])
if'LON_NUM_COEFF'in d:
lf.lon_num =list(map(float, d['LON_NUM_COEFF'].split()))
lf.lon_den =list(map(float, d['LON_DEN_COEFF'].split()))
lf.lat_num =list(map(float, d['LAT_NUM_COEFF'].split()))
lf.lat_den =list(map(float, d['LAT_DEN_COEFF'].split()))
elif dict_format =="rpcm":
lf.__dict__ = d
el:
rai ValueError(
"dict_format '{}' not supported. "
"Should be {{'geotiff','rpcm'}}".format(dict_format)
)
def projection(lf, lon, lat, alt):
"""
Convert geographic coordinates of 3D points into image coordinates.正投影：从地理坐标到图像坐标
Args:
lon (float or list): longitude(s) of the input 3D point(s)
lat (float or list): latitude(s) of the input 3D point(s)
alt (float or list): altitude(s) of the input 3D point(s)
Returns:
float or list: horizontal image coordinate(s) (column index, ie x)
float or list: vertical image coordinate(s) (row index, ie y)
"""
nlon =(np.asarray(lon)- lf.lon_offt)/ lf.lon_scale
nlat =(np.asarray(lat)- lf.lat_offt)/ lf.lat_scale
nalt =(np.asarray(alt)- lf.alt_offt)/ lf.alt_scale
col = apply_l_num, lf.col_den, nlat, nlon, nalt)
row = apply_w_num, lf.row_den, nlat, nlon, nalt)
col = col * lf.col_scale + lf.col_offt
row = row * lf.row_scale + lf.row_offt
return col, row
def localization(lf, col, row, alt, return_normalized=Fal):
"""
Convert image coordinates plus altitude into geographic coordinates.反投影：从图像坐标到地理坐标
Args:
col (float or list): x image coordinate(s) of the input point(s)
row (float or list): y image coordinate(s) of the input point(s)
alt (float or list): altitude(s) of the input point(s)
Returns:
float or list: longitude(s)
float or list: latitude(s)
"""
ncol =(np.asarray(col)- lf.col_offt)/ lf.col_scale
nrow =(np.asarray(row)- lf.row_offt)/ lf.row_scale
nalt =(np.asarray(alt)- lf.alt_offt)/ lf.alt_scale
if not hasattr(lf,'lat_num'):
lon, lat = lf.localization_iterative(ncol, nrow, nalt)
el:
lon = apply_rfm(lf.lon_num, lf.lon_den, nrow, ncol, nalt)
lat = apply_rfm(lf.lat_num, lf.lat_den, nrow, ncol, nalt)
流行的英文名if not return_normalized:
lon = lon * lf.lon_scale + lf.lon_offt
lat = lat * lf.lat_scale + lf.lat_offt
return lon, lat
return lon, lat
def localization_iterative(lf, col, row, alt):
"""
Iterative estimation of the localization function (image to ground),
for a list of image points expresd in image coordinates.
逆投影时的迭代函数
Args:
col, row: normalized image coordinates (between -1 and 1)
alt: normalized altitude (between -1 and 1) of the corresponding 3D
point
Returns:
lon, lat: normalized longitude and latitude
英语六级真题下载oyamaRais:
MaxLocalizationIterationsError: if the while loop exceeds the max
number of iterations, which is t to 100.
"""
# target point: Xf (f for final)
Xf = np.vstack([col, row]).T
# u 3 corners of the lon, lat domain and project them into the image
小学英语小报# to get the first estimation of (lon, lat)
# EPS is 2 for the first iteration, then 0.1.
lon =-col **0# vector of ones
lat =-col **0
EPS =2
x0 = apply_l_num, lf.col_den, lat, lon, alt)
y0 = apply_w_num, lf.row_den, lat, lon, alt)
x1 = apply_l_num, lf.col_den, lat, lon + EPS, alt)
y1 = apply_w_num, lf.row_den, lat, lon + EPS, alt)
x2 = apply_l_num, lf.col_den, lat + EPS, lon, alt)
y2 = apply_w_num, lf.row_den, lat + EPS, lon, alt)
n =0
while not np.all((x0 - col)**2+(y0 - row)**2<1e-18):
if n >100:
rai MaxLocalizationIterationsError("Max localization iterations (100) exceeded")
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
# project u on the ba (e1, e2): u = a1*e1 + a2*e2
# the exact computation is given by:
# M = np.vstack((e1, e2)).T
# a = np.dot(np.linalg.inv(M), u)
# but I don't know how to vectorize this.
# Assuming that e1 and e2 are orthogonal, a1 is given by
# <u, e1> / <e1, e1>
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()
# u the coefficients a1, a2 to compute an approximation of the
# point on the gound which in turn will give us the new X0
lon += a1 * EPS
lat += a2 * EPS
