我正在研究地球磁场方面的地面观测数据。我在Jupyter Notebook中使用Python。每个天文台都有一个数据序列,绘制成plt.scatter
图,如下所示:
验证码:
# Function that plots the time series and the data series from a .mat data file
# (in the dB-direction based on the out-variable), that is associated with the observatory code (based on the sta-variable)
# It also takes an optional argument 'showdatapoints' so its possible to show data points mainly for debugging
def plotdata_stations(filename, sta, out, start, stop, *args):
# Makes the function use one of the three columns (and thus directions of field) in the data series based on input
if out == 'radial':
dat = 0
elif out == 'theta':
dat = 1
elif out == 'phi':
dat = 2
# If the user inputs something else than one of these three, it returns an error and a hint to resolve the error
else:
print('\nError: Component ' + "'" + out + "'" + ' not recognized...')
print('\nHint: Try using ' + "'" + 'radial' + "', " + "'" + 'theta' + "' or " + "'" + 'phi' + "'.")
# Try to load file in case the file does not exist
try:
nam = DataLoad(filename, 'obs_all')
# If the file does not exist, it returns an error and a hint to resolve the error
except:
print('\nError: File ' + "'" + filename + "'" + ' not recognized...')
print('\nHint: Try putting your file in the same folder as this script.')
# If the observatory code is not in the data file, it returns an error and a hint to resolve the error
if sta not in nam and sta != 'all':
print('\nError: Observatory Code ' + "'" + sta + "'" + ' not recognized...')
print('\nHint: Try using an Observatory Code from the Observatory Locations map.')
# Load data from specific station and direction of field
dat = dB_stations(filename, sta, dat)
tim = t_stations(filename, sta)
# If start is a string but not 'min' and stop is a string but not 'max', print error
if (isinstance(start, str) == True and start != 'min') or (isinstance(stop, str) == True and stop != 'max'):
print('\nError: Time range is out of bounds in regards to data set...')
print('\nHint: Try using a time range that is included in data set or use ' + "'" + 'max' + "'" + ' or ' + "'" + 'min' + "'.")
# Set start and stop to min or max of time series of min or max is chosen
if start == 'min':
start = min(tim)
if stop == 'max':
stop = max(tim)
# Print error if time range is out of bounds in regards to the data set
if start < min(tim) or stop > max(tim):
print('\nError: Time range is out of bounds in regards to data set...')
print('\nHint: Try using a time range that is included in data set or use ' + "'" + 'max' + "'" + ' or ' + "'" + 'min' + "'.")
# Get all data associated with station in the specific time range
ndat = []
ntim = []
for x in range(0,len(tim)):
if start <= tim[x] <= stop:
ndat.append(dat[x])
ntim.append(tim[x])
# Sets up the data plot and shows it
plt.figure(figsize=(9,5))
plt.scatter(ntim, ndat, marker='.')
plt.title('Observatory Data from ' + sta + ' in the dB_' + out + '-direction')
plt.xlabel('Years')
plt.ylabel('nT/yr')
plt.show()
# For optional arguments
for x in args:
# If the optional argument 'showdatapoints' is called it will print the data points bellow the plot
if x == 'showdatapoints':
print('\nData series [nT/yr] in the dB_' + out + '-direction = \n')
print(dB_stations(filename, sta, dat))
print('\nTime series [Years] = \n')
print(t_stations(filename, sta))
# If the user inputs something else than 'showdatapoints', it returns an error and a hint to resolve the error
else:
print('\nError: Optional argument ' + "'" + x + "'" + ' not recognized...')
print('\nHint: Try using ' + "'" + 'showdatapoints' + "' or deleting the argument.")
我用“Cartopy”来绘制天文台的位置,就像这样:
验证码:
# Function that takes the name of a station, the direction of field, and the range of time to plot from a .mat data file
# and plots the data from the station at the location of the station
def staGlobal(filename, sta, out, start, stop):
# Load data
the = DataLoad(filename, 'theta_obs_all')
phi = DataLoad(filename, 'phi_obs_all')
nam = DataLoad(filename, 'obs_all')
# Uses removeRepetition from earlier on all three data sets
# Remember to use array[#,0] instead of just array[#] on theta and phi,
# as the theta and phi series is an array inside another array for some reason
the = removeRepetition(nam, the[:,0])
phi = removeRepetition(nam, phi[:,0])
nam = removeRepetition(nam, nam)
# Converts theta and phi from radians to degrees, as cartopy likes
the = np.rad2deg(the)
phi = np.rad2deg(phi)
# Converts theta from colatitude to latitude, as cartopy likes
the = -the+90
# Converts phi from range 0-360 to range -180-180, for convenience
for x in range(0,len(phi)):
if phi[x] >= 180:
phi[x] = phi[x]-360
# Load data from specific station and direction of field
dat = dB_stations(filename, sta, out)
tim = t_stations(filename, sta)
# Get location of station
for x in range(0,len(nam)):
if nam[x] == sta:
the = the[x]
phi = phi[x]
# Get all data associated with station in the specific time range
ndat = []
ntim = []
for x in range(0,len(tim)):
if start <= tim[x] <= stop:
ndat.append(dat[x])
ntim.append(tim[x])
# Trying to scale data down to point of obs (this part needs to be fixed/replaced)
###########################################
scax = 3
scay = 5
lonrat = scax*(25+max(tim)-max(ntim))
latrat = abs(min(ndat))
#if phi < 0:
# latrat = abs(min(ndat))-scay
#else:
# latrat = 0
ntim=np.multiply(ntim,scax)
ndat=np.divide(ndat,scay)
###########################################
# Plots data from obs at the location of obs
plt.scatter(phi, the, color='red', marker='.')
plt.scatter(ntim+phi+lonrat, ndat+the+latrat, color='blue', marker='.')
我所要做的就是在Map上的天文台的位置绘制出天文台的plt.scatter
图。我想让它看起来像这样:
是否有任何方法在Map上的特定点绘制散点图/曲线?如果你能帮我的话你绝对会保存我的项目。如果你需要更多的信息来帮助我,请随时询问。我是新来的这个论坛,所以请原谅我,如果这是太模糊。
1条答案
按热度按时间gev0vcfq1#
我是用数学方法做的:
但如果有人知道一个更“科学”的方法,我将欢迎你的答复。