pyrtklib -- 基于rtklib的python binding

香港理工大学IPNL实验室发布了基于rtklib的python binding -- pyrtklib(https://github.com/IPNL-POLYU/pyrtklib)。现在，你可以在python中自由调用rtklib头文件中所定义的类型与函数，由浅入深的使用rtklib的后处理，单点定位，RTK定位等功能。同时，在python中有大量易于使用的机器学习算法（如决策树,SVM,神经网络等），因此现在可以轻易的将机器学习算法与定位算法进行一些紧耦合的研究。python编写便于调试与敏捷开发，在算法测试完毕后，也基本可以1:1由C/C++重新复刻提高效率。在ros中也可配合rospy等编写收发节点。

这里是一份简短的中文教程。（本教程同样适用于基于rtklib_demo5的pyrtklib_demo5,只需要将pyrtklib替换为pyrtklib5)

1.安装

pip安装(推荐)

在linux与macos(intel)平台上，pyrtklib已经编译并上传至了pypi，因此，你可以直接使用pip来进行安装。

pip install pyrtklib

由于现在手上没有M芯片的macos，因此暂时没有办法编译适用于M芯片的版本。

手动安装

当然你也可以由源码安装

git clone https://github.com/IPNL-POLYU/pyrtklib

cd pyrtklib

python setup.py install # 可能需要sudo

注意，需要依赖cmake 3.0。没有在msvc上测试过，可能需要一些其他的设置。

2.基本特性

(a)使用rtklib中的结构体

rtklib中的所有结构体都已经被封装为了python类，因此可以通过类的方式来进行构造，比如：

import pyrtklib as prl 
obs = prl.obs_t() 
nav = prl.nav_t() 
sta = prl.sta_t()

这里的obs即是用于储存GNSS观测量的结构。所有类的成员也都被绑定，可以在python中直接进行取值与赋值操作（包括指针成员）。仍然以obs举例：

obs = prl.obs_t() 
obs.n = 10 
obs.nmax = 10 
obs.data = prl.Arr1Dobsd_t(10) #在下面会讲到数组

在obs结构体中，n与nmax代表这个观测文件中的储存数量与允许的最大数量，data字段是由obsd_t数组的指针，用于储存具体的伪距载波等数据。

(b)数组与指针

由于rtklib由C语言编写，大量使用函数副作用通过指针来传递返回值，因此需要在python中实现对其指针的封装。这也是这个binding的核心。这里采用了如下的结构来封装数组同时可以模拟指针：

template<class T>
struct Arr1D<T>{ 
    T* src; int len;
}

同时对它进行python绑定:

template<typename Type>
void bindArr1D(py::module_& m, const std::string& typeName) {
    auto cls = py::class_<Arr1D<Type>>(m, ("Arr1D" + typeName).c_str())
    .def(py::init([](int len){return std::unique_ptr<Arr1D<Type>>(new Arr1D<Type>(len));}))
    .def(py::init([](Type* src,int len){return std::unique_ptr<Arr1D<Type>>(new Arr1D<Type>((void*)src,len));}))
    .def("__len__",[](Arr1D<Type> &arr){return &arr.len;})
    .def("__getitem__",[](Arr1D<Type> &arr,int i){return arr.src+i;},py::return_value_policy::reference)
    .def("__getitem__",[](Arr1D<Type> &arr,py::slice slice){long start,stop,step;PySlice_Unpack(slice.ptr(),&start,&stop,&step);return new Arr1D<Type>(arr.src+start,stop-start);},py::return_value_policy::reference)
    .def("__setitem__",[](Arr1D<Type> &arr,int i, Type v){arr.src[i]=v;})
    .def("__iter__",[](Arr1D<Type> &arr){return pybind11::make_iterator(arr.src,arr.src+arr.len);})
    .def("deepcopy",static_cast<Arr1D<Type>*(Arr1D<Type>::*)()>(&Arr1D<Type>::deepcopy))
    .def("deepcopy",static_cast<Arr1D<Type>*(Arr1D<Type>::*)(int)>(&Arr1D<Type>::deepcopy))
    .def_readonly("ptr",&Arr1D<Type>::src,py::return_value_policy::reference)
    .def("set",[](Arr1D<Type> &arr,Arr1D<Type> *nsrc){arr.src = (Type*)nsrc->src;})
    .def("print",[](Arr1D<Type> &arr){std::cout<<(arr.src)<<std::endl;});

    if constexpr (std::is_same<Type, char>::value) {
        cls.def(py::init([](const std::string& s) {
                auto* arr = new Arr1D<char>(s.size()+1);
                std::memcpy(arr->src, s.data(), s.size());
                arr->src[s.size()] = '\0';  // Null-terminate the string
                return arr;
            }), py::arg("s"), "Constructor from Python str");
    }

}

在这里实现了两个构造函数，取序号与切片，赋值等操作，以上面的obsd_t为例。

import pyrtklib as prl
obsds = Arr1Dobsd_t(10) #构造一个长度为10的obsd_t数组
obsds_slice = Arr1Dobsd_t(obsds,5) #只含obsds中前5个元素的引用
obsds2 = obsds[3:8] #取切片
obsd = obsds[0] #取第一个元素，类型是obsd_t
obsds_list = list(obsds) #将数组变为python列表

在以上的构造函数中都不能通过提供一个列表来赋初值的操作，因为类型太复杂，需要考虑的问题太多，所以在构造后需要手动来进行赋值，比如：

import pyrtklib as prl
values = [1.1,1.2,1.3]
val = prl.Arr1Dfloat(3)
for i in range(3):
    val[i] = values[i]

唯一的例外是为字符串提供了这样一个构造函数：

import pyrtklib as prl
string = prl.Arr1Dchar("pyrtklib")

---

---

但是，在取值时请慎用！！！！！！这里并没有进行边界检查，因为有一些空间是由rtklib动态分配的，因此无法进行提前的边界检查，如果指针越界，将会造成段错误等后果！！！！！

---

3.简单应用

如果你只需要使用rtklib来获取最后的结果，那么可以直接使用postpos等方法，如下：

from pyrtklib import *

if __name__ == "__main__":
    gnss_file_list = []
    gnss_file_list.append("data/20210714.2.whampoa.ublox.f9p.obs")
    gnss_file_list.append("data/hksc1950.21f")
    gnss_file_list.append("data/hksc1950.21n")
    gnss_file_list.append("data/hksc1950.21o")
    output_path = Arr1Dchar("pyexample_output.txt")

    trace_lv = 3
    stat_lv = 3

    prcopt = prcopt_default
    solopt = solopt_default
    filopt = filopt_t()
    ts = gtime_t()
    te = gtime_t()
    tint = 0.0
    es = [2000,1,1,0,0,0]
    ee = [2000,12,31,23,59,59]
    pos = [0,0,0]
    i = 0
    j = 0
    n = 0
    ret = 0
    rov = ""
    base = ""
    ts.time = 0
    te.time = 0

    n = len(gnss_file_list)
    for i in gnss_file_list:
        print(i)

    print(output_path.ptr)

    snrmask = snrmask_t()
    snrmask.ena[0] = 1
    snrmask.ena[1] = 1
    for i in range(NFREQ):
        for j in range(9):
            snrmask.mask[i,j] = 10.0

    prcopt.mode = PMODE_KINEMA
    prcopt.soltype = 2
    prcopt.nf = 3
    prcopt.navsys = SYS_ALL

    prcopt.sateph = EPHOPT_BRDC
    prcopt.modear = 3
    prcopt.glomodear = 1
    prcopt.bdsmodear = 1
    prcopt.ionoopt = IONOOPT_BRDC
    prcopt.tropopt = TROPOPT_SAAS

    prcopt.posopt[4] = 0
    prcopt.rb[0] = -2414266.9197
    prcopt.rb[1] = 5386768.9868
    prcopt.rb[2] = 2407460.0314

    solopt.posf = SOLF_LLH
    solopt.times = TIMES_GPST
    solopt.timef = 0
    solopt.timeu = 3
    solopt.degf = 0
    solopt.outhead = 1
    solopt.outopt = 1

    solopt.datum = 0
    solopt.height = 0
    solopt.geoid = 0

    solopt.sstat = stat_lv
    solopt.trace = trace_lv

    print('starting rtklib ...')
    ret = postpos(ts,te,tint,0.0,prcopt,solopt,filopt,gnss_file_list,n,output_path,rov,base)
    print("rtklib status: %d"%ret)

这里在之前读入了接收机与基站的数据，然后进行了一些设置，最后使用postpos来进行结算，结果保存到了pyexample_output.txt。

4.深入应用

你可能需要深入的使用rtklib，这里给出一个单点定位的例子，我们分割出obs中的每一个历元片段，然后只对这个历元的单点定位进行求解。

from pyrtklib import *
import pandas as pd
import numpy as np


def nextobsf(obs,i):
    n = 0
    while i+n < obs.n:
        tt = timediff(obs.data[i+n].time,obs.data[i].time)
        if tt > 0.05:
            break
        n+=1
    return n

def get_sat_cnt(obs,i,m,rcv):
    cnt = 0
    for k in range(i,i+m):
        if obs.data[k].rcv == rcv:
            cnt+=1
    return cnt

def get_sat_pos(obsd,n,nav):
    svh = Arr1Dint(MAXOBS)
    rs = Arr1Ddouble(6*n)
    dts = Arr1Ddouble(2*n)
    var = Arr1Ddouble(1*n)
    satposs(obsd[0].time,obsd.ptr,n,nav,0,rs,dts,var,svh)
    noeph = []
    for i in range(n):
        if rs[6*i] == 0:
            noeph.append(i)
    nrs = Arr1Ddouble(6*(n-len(noeph)))
    j = 0
    for i in range(6*n):
        if rs[i]!=0:
            nrs[j] = rs[i]
            j+=1
    return nrs,noeph

def split_obs(obs):
    i = 0
    m = nextobsf(obs,i)
    obss = []
    while m!=0:
        tmp_obs = obs_t()
        tmp_obs.data = obs.data[i:i+m]
        tmp_obs.n = m
        tmp_obs.nmax = m
        i+=m
        obss.append(tmp_obs)
        m = nextobsf(obs,i)
    return obss
    
def get_obs_pnt(obs,nav,prcopt):
    m = obs.n
    sol = sol_t()
    sat = Arr1Dssat_t(MAXSAT)
    sol.time = obs.data[0].time
    msg = Arr1Dchar(100)
    azel = Arr1Ddouble(m*2)
    pntpos(obs.data.ptr,obs.n,nav,prcopt,sol,azel,sat.ptr,msg)
    return sol,sat,azel,msg

traceopen('test.log')
tracelevel(4)



files = [
    'data/KLB/F9P_211203_061807.obs',
    "data/KLB/BRDC00IGS_R_20213370000_01D_MN.rnx"
]
obs = obs_t()
nav = nav_t()
sta = sta_t()
prcopt = prcopt_default
solopt = solopt_default
readrnx(files[0],1,"",obs,nav,sta)
for i in files[1:]:
    readrnx(i,1,"",obs,nav,sta)
sortobs(obs)
uniqnav(nav)


data = []
obss = split_obs(obs)
prcopt.mode = PMODE_KINEMA
prcopt.navsys = SYS_ALL
prcopt.soltype = 0
prcopt.elmin = 0#15.0*D2R
prcopt.tidecorr = 0
prcopt.posopt[4] = 0
prcopt.tropopt = TROPOPT_SBAS
prcopt.ionoopt = IONOOPT_BRDC
prcopt.sateph = EPHOPT_BRDC
prcopt.modear = 3 #fix and hold

index = 0
for o in obss:
    sol,sat,azel,msg = get_obs_pnt(o,nav,prcopt)
    if msg[0] == "no observation data":
        break
    ecef = Arr1Ddouble(3)
    ecef[0] = sol.rr[0]
    ecef[1] = sol.rr[1]
    ecef[2] = sol.rr[2]
    pos = Arr1Ddouble(3)
    ecef2pos(ecef,pos)
    print(sol.time.time,*(np.array(pos[:2])*R2D),pos[2])
print('done')

这里我们调用了rtklib的readrnx，satposs等函数来读取与求解卫星位置等，最后使用pntpos只处理我们分割好的obs块。当然，你甚至可以更深入的自己重写pntpos，因为pntpos中的函数也全是由rtklib.h中的函数编写完成。

5.参考文献

如果您觉得这个工具有用，请引用我们的文章，万分感谢:)

@inproceedings{hu2023fisheye, title={Fisheye camera aided NLOS exclusion and learning-based pseudorange correction}, author={Hu, Runzhi and Wen, Weisong and Hsu, Li-ta}, booktitle={2023 IEEE 26th International Conference on Intelligent Transportation Systems (ITSC)}, pages={1--8}, year={2023}, organization={IEEE} }