您好,登录后才能下订单哦!
密码登录
登录注册
点击 登录注册 即表示同意《亿速云用户服务条款》
本篇内容主要讲解“Python如何根据站点列表绘制站坐标全球分布图”,感兴趣的朋友不妨来看看。本文介绍的方法操作简单快捷,实用性强。下面就让小编来带大家学习“Python如何根据站点列表绘制站坐标全球分布图”吧!
根据站点列表绘制站坐标全球分布图
输入:站点列表文件、SNX全球站点坐标文件
站点列表文件示例(可手动创建):
结果输出:
代码:
# 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如何根据站点列表绘制站坐标全球分布图”有了更深的了解,不妨来实际操作一番吧!这里是亿速云网站,更多相关内容可以进入相关频道进行查询,关注我们,继续学习!
免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:is@yisu.com进行举报,并提供相关证据,一经查实,将立刻删除涉嫌侵权内容。