如何用MLT替代地理纬度绘制南极上空MLAT-MLT能量分布极图？

问题核心是把地理坐标的海岸线转换为MLAT-MLT坐标，再在适配的坐标网格上绘图，而非依赖地理投影的worldmap。

关键步骤 #

1. 安装地磁坐标转换工具 #

MATLAB无内置MLAT-MLT转换，需借助第三方工具包，比如从MathWorks File Exchange下载的igrf工具包，或Space Physics Toolbox (SPT)。以SPT为例，先添加路径：

addpath('path/to/SPT');

2. 转换海岸线至MLAT-MLT坐标 #

MLT依赖观测时间，需与你的数据时间匹配，假设时间为obs_time = datetime('2023-01-01 12:00:00')：

load coastlines
% 筛选南极区域海岸线（纬度<-40）
south_idx = coastlat < -40;
clat = coastlat(south_idx);
clon = coastlon(south_idx);

% 地理坐标转地磁坐标
[mlat_coast, mlt_coast] = geo2mag(clat, clon, obs_time);

3. 绘制MLAT-MLT极坐标图 #

all_data = load('struct.mat');
all_data = all_data.struct;
lat = all_data.lat;
mlt = all_data.mlt;
heat = all_data.heat;

figure('Position', [100 100 800 800])
ax = polaraxes;
% MLT转弧度，MLAT转半径（-90对应0，-40对应50）
theta = mlt * 2*pi/24;
r = 90 + lat;
scatter(theta, r, 75, heat, 'filled', 's');

% 叠加转换后的海岸线
theta_coast = mlt_coast * 2*pi/24;
r_coast = 90 + mlat_coast;
plot(theta_coast, r_coast, 'k', 'LineWidth', 1);

% 美化设置
ax.ThetaDir = 'clockwise';  % MLT顺时针递增
ax.ThetaZeroLocation = 'top';  % 0MLAT置于顶部
ax.RTick = [10 30 50];
ax.RTickLabel = {'-80', '-60', '-40'};
ax.ThetaTick = 0:3:24;
ax.ThetaTickLabel = {'0', '3', '6', '9', '12', '15', '18', '21', '24'};
colormap('parula')
caxis([0 40000])
colorbar;
title('南极离子能量分布 (MLAT vs MLT)');

Python用pyspedas做地磁转换，matplotlib绘图，灵活性更强。

关键步骤 #

1. 安装依赖 #

pip install pyspedas matplotlib cartopy

2. 代码实现 #

import numpy as np
import matplotlib.pyplot as plt
from pyspedas import geomag
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader

# 加载数据
all_data = np.load('struct.mat', allow_pickle=True)['struct'].item()
lat = all_data['lat']
mlt = all_data['mlt']
heat = all_data['heat']

# 匹配数据的观测时间
obs_time = '2023-01-01 12:00:00'

# 转换南极海岸线至MLAT-MLT
land_feature = cfeature.NaturalEarthFeature('physical', 'land', '110m',
                                            edgecolor='black', facecolor='none')
reader = Reader(land_feature.path)
mlat_coast = []
mlt_coast = []

for geom in reader.geometries():
    for polygon in geom.geoms:
        coords = np.array(polygon.exterior.coords)
        clon, clat = coords[:, 0], coords[:, 1]
        south_idx = clat < -40
        if np.any(south_idx):
            # 地理坐标转地磁坐标
            mlat, _, mlt_val = geomag.geomag(clat[south_idx], clon[south_idx], obs_time, 'geo', 'mag')
            mlat_coast.extend(mlat)
            mlt_coast.extend(mlt_val)

# 转换为极坐标格式
theta = np.deg2rad(mlt * 15)
r = 90 + lat
theta_coast = np.deg2rad(np.array(mlt_coast) * 15)
r_coast = 90 + np.array(mlat_coast)

# 绘图
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, polar=True)
sc = ax.scatter(theta, r, c=heat, s=75, marker='s', cmap='viridis', vmin=0, vmax=40000)
ax.plot(theta_coast, r_coast, 'k', linewidth=1)

# 美化设置
ax.set_theta_direction('clockwise')
ax.set_theta_zero_location('N')
ax.set_xticks(np.deg2rad(np.arange(0, 24, 3)*15))
ax.set_xticklabels([str(i) for i in range(0,24,3)] + ['24'])
ax.set_yticks([10, 30, 50])
ax.set_yticklabels(['-80', '-60', '-40'])
plt.colorbar(sc, label='Energy')
plt.title('南极离子能量分布 (MLAT vs MLT)')
plt.show()

注意事项 #

• 地磁转换必须指定准确的观测时间，MLT是时间依赖参数，不同时间海岸线在MLAT-MLT图上的位置会旋转，与参考图效果一致。
• MATLAB若无法获取SPT，可使用igrf12syn函数手动计算地磁坐标。

内容的提问来源于stack exchange，提问作者steve