Python根据站点列表绘制站坐标全球分布图的示例

根据站点列表绘制站坐标全球分布图
输入：站点列表文件、snx全球站点坐标文件
站点列表文件示例(可手动创建)：

snx全球站点坐标文件下载地址：
ftp://igs.gnsswhu.cn/pub/whu/pub/gps/products/yyyy/igsyypwwww.snx.z
结果输出：

代码:

# coding=utf-8
# !/usr/bin/env python
'''
 program:plot_global_sitemap.py 
 function:根据站点列表绘制站坐标全球分布图
 author:lz_cumt
 version:1.0
 date:2021/12/10
 '''
from math import pi, sqrt, atan, atan2, sin, cos
import matplotlib.pyplot as plt
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import longitudeformatter, latitudeformatter

# xyz转换为llh（经纬度）
def xyz2llh(ecef, site):
    aell = 6378137.0
    fell = 1.0 / 298.257223563
    deg = pi / 180
    u = ecef[0]
    v = ecef[1]
    w = ecef[2]
    esq = 2*fell-fell*fell
    lat = 0
    n = 0
    if w == 0:
        lat = 0
    else:
        lat0 = atan(w/(1-esq)*sqrt(u*u+v*v))
        j = 0
        delta = 10 ^ 6
        limit = 0.000001/3600*deg
        while delta > limit:
            n = aell / sqrt(1 - esq * sin(lat0)*sin(lat0))
            lat = atan((w / sqrt(u*u + v*v)) * (1 + (esq * n * sin(lat0) / w)))
            delta = abs(lat0 - lat)
            lat0 = lat
            j = j + 1
            if j > 10:
                break
    long = atan2(v, u)
    h = (sqrt(u*u+v*v)/cos(lat))-n
    llh = [site, long * 180 / pi, lat * 180 / pi, h]
    return llh

# 由站点文件获取站点列表存入sitelist
def getsite(listfile):
    sitelist = []
    f = open(listfile)
    ln = f.readline()
    while ln:
        sitelist.append(ln[0:4].upper())
        ln = f.readline()
    return sitelist

# 根据站点名在snx文件中搜索xyz坐标转化为经纬度并输出
def getblh_single(site,snxlines):
    xyz = [0, 0, 0]
    for ln in snxlines:
        if site in ln:
            if 'stax   ' in ln:
                xyz[0] = float(ln[47:68])
            if 'stay   ' in ln:
                xyz[1] = float(ln[47:68])
            if 'staz   ' in ln:
                xyz[2] = float(ln[47:68])
    blh = xyz2llh(xyz, site)
    if len(blh) != 4:
        print('[info] sitecrd for', site, 'is not found in the snxfile')
    return blh

def getblh(listfile, snxfile):
    siteblh = []
    sitelist = getsite(listfile)
    f = open(snxfile)
    lns = f.readlines()
    for site in sitelist:
        siteblh.append(getblh_single(site, lns))
    return siteblh

def plotsite(siteblh):
    # mpl.rcparams['font.sans-serif'] = ['helvetical']
    mpl.rcparams['axes.unicode_minus'] = false
    mpl.rc('xtick', labelsize=9)
    mpl.rc('ytick', labelsize=9)
    mpl.rcparams['xtick.direction'] = 'in'
    mpl.rcparams['ytick.direction'] = 'in'

    fig = plt.figure(figsize=(14, 7))
    ax = plt.axes(projection=ccrs.platecarree(central_longitude=150))
    ax.set_extent([-180, 180, -90, 90], crs=ccrs.platecarree())
    ax.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.platecarree())
    ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.platecarree())
    ax.add_feature(cfeature.land)
    ax.add_feature(cfeature.ocean)
    ax.add_feature(cfeature.coastline, linewidth=0.1)

    for site in siteblh:
        ax.plot(site[1], site[2], 'o', color='r', mec='k', mew=0.5, transform=ccrs.geodetic(), ms=13.0)
        plt.text(site[1] + 1.5, site[2] + 1.5, site[0], transform=ccrs.geodetic(),fontsize='x-large')  # 添加站名标注
    plt.xticks(fontsize='x-large')
    plt.yticks(fontsize='x-large')
    lon_formatter = longitudeformatter(zero_direction_label=true)
    lat_formatter = latitudeformatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)

    fig.savefig('global_sitemap.png', bbox_inches='tight', dpi=400)
    plt.show()


if __name__ == '__main__':
    listfile = r'site.info'               # 输入要画的站点列表文件
    snxfile = r'igs21p2177.snx'              # 输入igs站坐标文件
    siteblh = getblh(listfile, snxfile)   # 获取所有站点的经纬度
    plotsite(siteblh)           # 画图
    print('[info] plot complete!')        # 完成

到此这篇关于python根据站点列表绘制站坐标全球分布图的文章就介绍到这了,更多相关python绘制站坐标全球分布图内容请搜索以前的文章或继续浏览下面的相关文章希望大家以后多多支持！