如何在cartopy中绘制地图上特定点的散点图

Question

我正在研究有关地球磁场的地面观测数据。我在 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='.')

我想做的就是在地图上的天文台位置绘制天文台的

plt.scatter

图。我想让它看起来像这样：

有没有办法在地图上的特定点绘制散点图/曲线？ 如果你能帮助我，你绝对会拯救我的项目。如果您需要更多信息来帮助我解决这个问题，请随时询问。我是这个论坛的新手，如果这太含糊，请多多包涵。

Answer 1

我是这样用数学方法做的：

# Scaling data down to point of obs
###########################################
scax = 3
scay = 4
ntim=np.subtract(ntim,2000)
ntim=np.multiply(ntim,scax)
ndat=np.divide(ndat,scay)
###########################################

# Plots data from obs at the location of obs
# Plot at locations by adding coordinates of obs and subtracting the mean of each axis of data set
plt.scatter(ntim-np.mean(ntim)+phi, ndat-np.mean(ndat)+the, color='blue', marker='.')
plt.scatter(phi, the, color='red', marker='.')

但如果有人知道更“科学”的方法，我欢迎您的回复。