可视卫星数与PDOP图绘制

---

前言

GNSS相关的论文中，经常涉及可视卫星数以及PDOP值的绘制，虽然有时有水字数之嫌，但不少情况下也是必须的。这里基于rtklib的函数，写了一下方便自定义输出的计算函数。

---

一、主要思路和需要调用的函数

1、首先写一个读取obs文件的函数，只读取卫星数和卫星列表即可。需要的话，也可以读取概略坐标。
2、计算卫星数，简单利用satsys判断一下就行了。
因为我需要考虑PPP-B2b匹配的问题，因此作了更多的判断。
3、计算pdop，这个需要调用的函数略多一点，封装了一个新的函数
主要涉及的rtklib原函数如下：

readrnxnav()	//读取星历
ephpos()		//计算卫星坐标
geodist()		//这个只是为了辅助计算方位角
satazel()		//方位角
dops()			//用方位角算dops

二、代码实现和相关的函数

int main(int argc, char** argv)
{
	//这里实现自动化会更好一点
	char pathObs[100] = "D:\\D-PaperData\\1210-orbit\\cxws3440.23o";
	char pathLNAV[100] = "D:\\D-PaperData\\1210-orbit\\BRDC00IGS_R_20233440000_01D_MN.rnx";
	//下面两行涉及PPP-B2b的问题，如果不涉及，直接删去就好。后面B2b相关的也是一样，不再赘述
	char pathCNAV[100] = "D:\\D-PaperData\\1210-orbit\\tarc3440.23b_cnav";
	char pathB2b[100] = "D:\\D-PaperData\\fixfiles\\PPPB2b";
	char pathOut[100] = "D:\\D-PaperData\\1210-orbit\\sumPdop0121.txt";

	gtime_t ts = { 0 }, te = { 0 }, tmid = { 0 };
	int i, j, k, nums[4], BnG, BnC, stat1, svh[60], satG[32], satC[60];
	double tss[6] = { 2023,12,10,2,0,0 };
	double tee[6] = { 2023,12,10,14,0,0 };
	ts = epoch2time(tss);
	te = epoch2time(tee);
	double ti = 6.0, pG, pC,pGC;
	navs.eph = NULL; navs.n = navs.nmax = 0;
	char outbuff[200], strtime[50];
	double rs[6 * 100], dts[2 * 100], var[100], azel[2 * 100], rb[3];
	int bds3 = satid2no("C19"),ind;

	//读取数据
	ReadPPPB2b(pathB2b, ts, te);
	new_readrnxnav(pathLNAV, &navs, "-SYS=G");
	//readbrd4(pathCNAV, &navs, "-SYS=C");
	new_readrnxnav(pathCNAV, &navs, "-SYS=C");
	ReadObsSat(pathObs, nsats, rb,ts,te); 
	std::cout << "read NAV number: " << navs.n << "\n";
	std::cout << "read OBS number: " << nsats.size() << "\n";

	//打开和输出文件头
	FILE* fpout;
	if (!(fpout = fopen(pathOut, "w"))) {
		return 0;
	}
	fprintf(fpout, "       time         nG nC nE nR BnG BnC pdopG pdopC pdopGC\n");   //b表示PPPB2b，p即PDOP

	//循环处理数据
	for (i = 0; i < nsats.size(); i++)
	{
		nobs_t ttobs = nsats[i];
		UpdateB2b(ttobs.t0, navs.b2bs); //更新B2b，不涉及可以直接删去，这里就不放原型了
		std::memset(nums, 0, sizeof(int) * 4); //初始化为0
		BnG = BnC = pG = pC = pGC = 0;
		//这里只考虑四个系统，统计卫星数
		for (j = 0; j < ttobs.n; j++) {
			switch (satsys(ttobs.satno[j], NULL)) {
			case SYS_GPS: satG[nums[0]] = ttobs.satno[j]; nums[0]++; break;   //存卫星号
			case SYS_CMP: if (ttobs.satno[j] >= bds3) { satC[nums[1]] = ttobs.satno[j]; nums[1]++; } break;
			case SYS_GAL: nums[2]++; break;
			case SYS_GLO: nums[3]++; break;
			default: break;
			}
		}
		//计算pdop和PPP-B2b ,原函数为ephpos，几乎不用改动
		for (j = 0; j < nums[0]; j++) {   //GPS
			stat1 = ephpos_b2b(ttobs.t0, ttobs.t0, satG[j], &navs, 1, rs + j * 6, dts + j * 2, var + j, svh + j);			
			if (stat1) BnG++;
		}
		pG = CalPdop(nums[0], rs, rb);
		for (; j < nums[0]+nums[1]; j++) {   //BDS
			ind = j - nums[0];
			stat1 = ephpos_b2b(ttobs.t0, ttobs.t0, satC[ind], &navs, 1, rs + j * 6, dts + j * 2, var + j, svh + j);			
			if (stat1) BnC++;
		}	
		pC = CalPdop(nums[1], rs + nums[0] * 6, rb);
		pGC = CalPdop(nums[0] + nums[1], rs, rb);
		time2str(ttobs.t0, strtime, 0);
		fprintf(fpout, "%s %2d %2d %2d %2d %2d  %2d  %4.1f  %4.1f  %4.1f\n", strtime, nums[0]
			, nums[1], nums[2], nums[3], BnG, BnC, pG, pC, pGC);

		sprintf(outbuff, "\rprocess: %d/%d", i, nsats.size());
		cout << outbuff;

	}
	fclose(fpout);
}

//下面是一些添加的函数
//存储当前历元的时间、历元、卫星号
typedef struct {
	int n;
	int satno[MAXSAT];
	gtime_t t0;
} nobs_t;
vector<nobs_t> nsats;  //一个元素就是一个观测历元

//利用观测文件读取obs数据，主要是卫星数和列表
static int ReadObsSat(const char* file, vector<nobs_t>&nsats, double *rb,gtime_t ts, gtime_t te)
{
	char buff[MAXREADLEN],satid[4];
	FILE* fp;
	gtime_t t0;
	int nobs,i,satno;
	if (!(fp = fopen(file, "r"))) {
		return 0;
	}
	while (fgets(buff, MAXREADLEN, fp)) {   //计算pdop，用概略位置就行
		if (strstr(buff, "APPROX POSITION XYZ")) {
			for (i = 0; i < 3; i++) rb[i] = str2num(buff, i * 14, 14);
		}
		if (strstr(buff, "END OF HEADER")) break;
	}	
	while (fgets(buff, MAXREADLEN, fp)) {
		if (buff[0] != '>') continue;
		nobs_t ttobs;
		str2time(buff, 1, 27, &t0);
		if (timediff(t0, ts) < 0) continue;
		if (timediff(t0, te) > 0) break;
		nobs=(int)str2num(buff, 33, 2);
		ttobs.n = nobs;
		ttobs.t0 = t0;
		for (i = 0; i < nobs; i++) {
			fgets(buff, MAXREADLEN, fp);
			strncpy(satid, buff, 3);
			satid[3] = '\0';
			satno=satid2no(satid);
			ttobs.satno[i] = satno;
		}
		nsats.push_back(ttobs);
	}
	return 1;
}

//计算pdop,卫星数,卫星坐标,和移动站坐标;和rtklib格式是一样的
static double CalPdop(const int num, double *rs, double *rb) {
	double pos[3],e[3],*azel,dop[4];	
	double minele = PI / 180 * 7;
	azel = mat(2, num);
	ecef2pos(rb, pos);
	for (int i=0; i < num; i++) {
		geodist(rs+i*6,rb,e);
		satazel(pos, e, azel+i*2);		
	}
	dops(num, azel, minele, dop);
	return dop[1];
}

三、画图

这里输出后用python画的图

import matplotlib.pyplot as plt
import numpy as np
import datetime
from proplot import rc
import matplotlib.dates as mdates

def ReadTxt(path):
    data=[]
    times=[]
    fileHandler=open(path,  "r")  
    line  =  fileHandler.readline() #去掉文件头
    while True:
        line  =  fileHandler.readline()
        if(not line) or line.find("EOF")>-1:
            break;
        s=line.split()
        t0=datetime.datetime.strptime(s[0]+' '+s[1], '%Y/%m/%d %H:%M:%S')
        nG=float(s[2])
        nC=float(s[3])
        BnG=float(s[6])
        BnC=float(s[7])
        pG=float(s[8])
        pC=float(s[9])
        pGC=float(s[10])
        data.append([nG,nC,BnG,BnC,pG,pC,pGC])
        times.append(t0)
    data=np.array(data)
    return times,data

if __name__=="__main__": 
    # In[0] 读取数据
    path='sumPdop0121.txt' 
    datestr='2024/01/21'   
    ts=datetime.datetime.strptime(datestr+' 02:00:00.0', '%Y/%m/%d %H:%M:%S.%f')
    te=datetime.datetime.strptime(datestr+' 10:00:00.0', '%Y/%m/%d %H:%M:%S.%f') 
    times,data=ReadTxt(path)
    
    # In[1]画图
    # 统一设置字体
    rc["font.family"] = "TeX Gyre Schola"
    # 统一设置轴刻度标签的字体大小
    rc['tick.labelsize'] = 13
    # 统一设置xy轴名称的字体大小
    rc["axes.labelsize"] = 13
    # 统一设置轴刻度标签的字体粗细
    rc["axes.labelweight"] = "light"
    rc['xtick.direction'] = 'in'
    rc['ytick.direction'] = 'in'
      
    fig, axs = plt.subplots(2, 1, sharex=True) 
    fig.set_facecolor('white')
    fig.set_size_inches(8, 4) 
    cs=['#2C906A','#C7590C','#1F64A5','#A01B2D']   #绿 蓝 橙 暗红
         
    axs[0].plot(times,data[:,0],label='GPS',color=cs[0])
    axs[0].plot(times,data[:,1],label='BDS',color=cs[1])
    axs[0].plot(times,data[:,2],label='PPP-B2b GPS',color=cs[2])
    axs[0].plot(times,data[:,3],label='PPP-B2b BDS',color=cs[3])
    axs[0].set_ylim((0, 17))  
    axs[0].set_ylabel('NSAT')
    axs[0].legend(loc='lower left',ncol=2,bbox_to_anchor=(0.05, 0.05))

    axs[1].plot(times,data[:,4],label='GPS',color=cs[0])
    axs[1].plot(times,data[:,5],label='BDS',color=cs[1])
    axs[1].plot(times,data[:,6],label='GPS+BDS',color=cs[2])
    axs[1].legend(loc='lower left',ncol=3,bbox_to_anchor=(0.05, 0.05))
    #设置x轴
    timeFmt = mdates.DateFormatter('%H:%M') 
    axs[1].xaxis.set_major_formatter(timeFmt)
    ts+=datetime.timedelta(seconds=-300)
    te+=datetime.timedelta(seconds=300)
    axs[1].set_xlim(ts, te)
    axs[1].set_ylim((0, 5.1))  
    axs[1].set_ylabel('PDOP')
    
    # 次要刻度 网格
    for row in range(2):  # 遍历每一行
        axs[row].tick_params(axis='x',direction='in',which='minor', bottom=False)  
        axs[row].tick_params(axis='y',direction='in',which='minor', left=False, right=False)
        axs[row].xaxis.set_minor_locator(mdates.HourLocator(interval=2))
        axs[row].grid(False)
    plt.subplots_adjust(left=0.08,right=0.98,top=0.95,bottom=0.10,hspace=0.1,wspace=0.05)    

四、示例

输出格式如下，根据自己的需要改改：

       time         nG nC nE nR BnG BnC pdopG pdopC pdopGC
2024/01/21 02:00:00  9 13  5  6  8   8   1.8   2.1   1.3
2024/01/21 02:00:01  9 13  5  6  8   8   1.8   2.1   1.3
2024/01/21 02:00:02  9 13  5  6  8   8   1.8   2.1   1.3
2024/01/21 02:00:03  9 13  5  6  8   8   1.8   2.1   1.3
2024/01/21 02:00:04  9 13  5  6  8   8   1.8   2.1   1.3

绘图结果如下，当然不是很完美：