How to draw the global distribution map of station coordinates according to the site list by Python
This article mainly explains "Python how to draw the global distribution map of station coordinates according to the site list", interested friends may wish to have a look. The method introduced in this paper is simple, fast and practical. Next let the editor to take you to learn "Python how to draw the global distribution map of station coordinates according to the site list"!
Draw the global distribution map of station coordinates according to the list of sites
Input: site list file, SNX global site coordinate file
Example of a site list file (which can be created manually):
Result output:
Code:
# coding=utf-8#! / usr/bin/env python''' Program:plot_global_sitemap.py Function: according to the list of sites, draw the global distribution map of station coordinates Author:LZ_CUMT Version:1.0 Date:2021/12/10''from math import pi, sqrt, atan, atan2, sin, cosimport matplotlib.pyplot as pltimport matplotlib as mplimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter# xyz into llh (latitude and longitude) def xyz2llh (ecef) Site): aell = 6378137.0 fell = 6378137.0 / 298.257223563 deg = pi / 180U = ecef [0] v = ecef [1] w = ecef [2] esq = 2*fell-fell*fell lat = 0N = 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 (ubiu + vav)) * (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# gets the site list from the site file and stores it in sitelistdef getSite (listfile): sitelist = [] f = open (listfile) ln = f.readline () while ln: sitelist.append (ln [0:4] .upper ()) ln = f.readline () return sitelist# searches the snx file for XYZ coordinates and converts them to longitude and latitude and outputs 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 blhdef getBLH (listfile, snxfile): siteBLH = [] sitelist = getSite (listfile) f = open (snxfile) lns = f.readlines () for site in sitelist: siteBLH.append (getBLH_single (site Lns) return siteBLHdef 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) ) 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,30,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], 'oasis, 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') # add the site name label 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 = ruplosite.info' # enter the site list file snxfile = ruploigs21P2177.snx' # enter the IGS station coordinate file siteBLH = getBLH (listfile Snxfile) # get the longitude and latitude of all sites plotsite (siteBLH) # drawing print ('[INFO] Plot completions') # this is done I believe you have a deeper understanding of "Python how to draw the global distribution map of station coordinates according to the site list". You might as well do it in practice. Here is the website, more related content can enter the relevant channels to inquire, follow us, continue to learn!
