/
cosmic_and_ivm_demo.py
354 lines (276 loc) · 11.9 KB
/
cosmic_and_ivm_demo.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
"""Seasonal analysis demo using COSMIC RO profiles and IVM in situ data.
"""
import datetime as dt
import matplotlib.pyplot as plt
import numpy as np
import numpy.ma as ma
from packaging import version as pack_version
import pandas as pds
from scipy.stats import mode
import apexpy
import pysat
import pysatCDAAC
import pysatNASA
import pysatSeasons
if pack_version(pysatCDAAC.__version__) <= pack_version('0.0.2'):
estr = ' '.join(['This demo is written expecting xarray support for ',
'COSMIC data. Unfortunately, this is not supported by ',
'the currently installed version. Please see the demo ',
'code in pysatSeasons v0.1.3 for COSMIC support when in ',
'pandas data format.'])
raise(ValueError, estr)
def add_magnetic_coordinates(inst):
"""Add magnetic longitude and latitude for COSMIC density maximum location.
Adds 'qd_lon' and 'qd_lat' to `inst`.
Parameters
----------
inst : pysat.Instrument
'COSMIC' Instrument object
"""
apex = apexpy.Apex(date=inst.date)
# Convert geographic profile location to magnetic location
lats, lons = apex.geo2qd(inst['edmaxlat'], inst['edmaxlon'],
inst['edmax_alt'])
# Longitudes between 0 - 360.
idx, = np.where(lons < 0)
lons[idx] += 360.
# Add data and metadata to Instrument object
inst['edmax_qd_lat'] = lats
notes_str = ''.join(['Obtained from apexpy by transforming ',
'`edmaxlat` and `edmaxlon` into quasi-dipole ',
'coordinates.'])
meta_data = {inst.meta.labels.units: 'degrees',
inst.meta.labels.name: 'Quasi-Dipole Latitude',
inst.meta.labels.notes: notes_str,
inst.meta.labels.fill_val: np.nan,
inst.meta.labels.min_val: -90.,
inst.meta.labels.max_val: 90.}
inst.meta['edmax_qd_lat'] = meta_data
inst['edmax_qd_lon'] = lons
meta_data = {inst.meta.labels.units: 'degrees',
inst.meta.labels.name: 'Quasi-Dipole Longitude',
inst.meta.labels.notes: notes_str,
inst.meta.labels.fill_val: np.nan,
inst.meta.labels.min_val: 0.,
inst.meta.labels.max_val: 360.}
inst.meta['edmax_qd_lon'] = meta_data
return
def restrict_abs_values(inst, label, max_val):
"""Restrict absolute quasi-dipole latitude values.
Parameters
----------
inst : pysat.Instrument
'COSMIC' Instrument object
label : str
Label for variable to restrict.
max_val : float
Absolute maximum value of `label`. Values greater
than `max_val` are removed from `inst`.
"""
inst.data = inst[np.abs(inst[label]) <= max_val]
return
def filter_values(inst, label, val_range=(0., 15.)):
"""Filter values to those in `label` that are within `val_range` limits.
Parameters
----------
inst : pysat.Instrument
'COSMIC' Instrument object
label : str
Label for variable to restrict.
val_range : tuple/list of floats
Minimum and maximum values allowed. Times where `label`
outside `val_range` are removed from `inst`. (default=(0., 15.))
"""
inst.data = inst[(inst[label].values >= val_range[0])
& (inst[label].values <= val_range[1])]
return
def add_log_density(inst):
"""Add the log of COSMIC maximum density, 'edmax'.
Parameters
----------
inst : pysat.Instrument
'COSMIC' 'GPS' object.
"""
# Assign data
inst['lognm'] = np.log10(inst['edmax'])
# Assign metadata
notes_str = ''.join(['Log base 10 of `edmax` variable.'])
meta_data = {inst.meta.labels.units: 'Log(N/cc)',
inst.meta.labels.name: 'Log Density',
inst.meta.labels.notes: notes_str,
inst.meta.labels.fill_val: np.nan,
inst.meta.labels.min_val: -np.inf,
inst.meta.labels.max_val: np.inf}
inst.meta['lognm'] = meta_data
return
def add_scale_height(inst):
"""Calculate topside scale height using observed electron density profiles.
Parameters
----------
inst : pysat.Instrument
'COSMIC' 'GPS' Instrument.
"""
output = inst['edmaxlon'].copy()
for i, profile in enumerate(inst['ELEC_dens']):
profile = profile[(profile
>= (1. / np.e) * inst[i, 'edmax'])
& (profile.coords["MSL_alt"] >= inst[i, 'edmax_alt'])]
# Want the first altitude where density drops below NmF2/e.
if len(profile) > 10:
i1 = profile.coords["MSL_alt"][1:]
i2 = profile.coords["MSL_alt"][0:-1]
mode_diff = mode(i1.values - i2.values)[0][0]
# Ensure there are no gaps, if so, remove all data above gap
idx, = np.where((i1.values - i2.values) > 2 * mode_diff)
if len(idx) > 0:
profile = profile[0:idx[0]]
# Ensure there are no null values in the middle of profile.
idx, = np.where(profile.isnull())
if len(idx) > 0:
profile = profile[0:idx[0]]
if len(profile) > 10:
# Make sure density at highest altitude is near Nm/e
if profile[-1] / profile[0] < 0.4:
alt = profile.coords["MSL_alt"]
alt_diff = alt[-1] - alt["MSL_alt"][0]
if alt_diff >= 500:
alt_diff = np.nan
output[i] = alt_diff
else:
output[i] = np.nan
else:
output[i] = np.nan
inst['thf2'] = output
return
# Register all instruments in pysatCDAAC and pysatNASA. Only required once
# per install.
pysat.utils.registry.register_by_module(pysatCDAAC.instruments)
pysat.utils.registry.register_by_module(pysatNASA.instruments)
# Dates for demo
ssn_days = 67
start_date = dt.datetime(2009, 12, 21) - pds.DateOffset(days=ssn_days)
stop_date = dt.datetime(2009, 12, 21) + pds.DateOffset(days=ssn_days)
# Instantiate IVM Object
ivm = pysat.Instrument(platform='cnofs', name='ivm', tag='',
clean_level='clean')
# Restrict measurements to those near geomagnetic equator.
ivm.custom_attach(restrict_abs_values, args=['mlat', 25.])
# Perform seasonal average
ivm.bounds = (start_date, stop_date)
ivm_results = pysatSeasons.avg.median2D(ivm, [0, 360, 24], 'alon',
[0, 24, 24], 'mlt',
['ionVelmeridional'])
# Create COSMIC instrument object. Engage supported keyword `altitude_bin`
# to bin all altitude profiles into 3 km increments.
cosmic = pysat.Instrument(platform='cosmic', name='gps', tag='ionprf',
clean_level='clean', altitude_bin=3)
# Apply custom functions to all data that is loaded through cosmic
cosmic.custom_attach(add_magnetic_coordinates)
# Select locations near the magnetic equator
cosmic.custom_attach(filter_values, args=['edmax_qd_lat', (-10., 10.)])
# Take the log of NmF2 and add to the dataframe
cosmic.custom_attach(add_log_density)
# Calculates the height above hmF2 to reach Ne < NmF2/e
cosmic.custom_attach(add_scale_height)
# Perform a bin average of multiple COSMIC data products, from start_date
# through stop_date. A mixture of 1D and 2D data is averaged.
cosmic.bounds = (start_date, stop_date)
cosmic_results = pysatSeasons.avg.median2D(cosmic, [0, 360, 24], 'edmax_qd_lon',
[0, 24, 24], 'edmaxlct',
['ELEC_dens', 'edmax_alt',
'lognm', 'thf2'])
# The work is done, plot the results!
# Make IVM and COSMIC plots
fig, axarr = plt.subplots(4, sharex=True, sharey=True, figsize=(8.5, 11))
cax = []
# Meridional ion drift average
mer_drifts = ivm_results['ionVelmeridional']['median']
x_arr = ivm_results['ionVelmeridional']['bin_x']
y_arr = ivm_results['ionVelmeridional']['bin_y']
# Mask out NaN values
masked = np.ma.array(mer_drifts, mask=np.isnan(mer_drifts))
# Plot, NaN values are white.
# Note how the data returned from the median function is in plot order.
cax.append(axarr[0].pcolor(x_arr, y_arr,
masked, vmax=30., vmin=-30.,
edgecolors='none'))
axarr[0].set_ylim(0, 24)
axarr[0].set_yticks([0, 6, 12, 18, 24])
axarr[0].set_xlim(0, 360)
axarr[0].set_xticks(np.arange(0, 420, 60))
axarr[0].set_ylabel('Magnetic Local Time (Hours)')
axarr[0].set_title('IVM Meridional Ion Drifts')
cbar0 = fig.colorbar(cax[0], ax=axarr[0])
cbar0.set_label('Ion Drift (m/s)')
max_dens = cosmic_results['lognm']['median']
cx_arr = cosmic_results['lognm']['bin_x']
cy_arr = cosmic_results['lognm']['bin_y']
# Mask out NaN values
masked = np.ma.array(max_dens, mask=np.isnan(max_dens))
# Plot, NaN values are white
cax.append(axarr[1].pcolor(cx_arr, cy_arr,
masked, vmax=6.1, vmin=4.8,
edgecolors='none'))
axarr[1].set_title('COSMIC Log Density Maximum')
axarr[1].set_ylabel('Solar Local Time (Hours)')
cbar1 = fig.colorbar(cax[1], ax=axarr[1])
cbar1.set_label('Log Density (N/cc)')
max_alt = cosmic_results['edmax_alt']['median']
# Mask out NaN values
masked = np.ma.array(max_alt, mask=np.isnan(max_alt))
# Plot, NaN values are white
cax.append(axarr[2].pcolor(cx_arr, cy_arr,
masked, vmax=375., vmin=200.,
edgecolors='none'))
axarr[2].set_title('COSMIC Altitude Density Maximum')
axarr[2].set_ylabel('Solar Local Time (Hours)')
cbar = fig.colorbar(cax[2], ax=axarr[2])
cbar.set_label('Altitude (km)')
max_th = cosmic_results['thf2']['median']
# Mask out NaN values
masked = np.ma.array(max_th, mask=np.isnan(max_th))
# Plot, NaN values are white
cax.append(axarr[3].pcolor(cx_arr, cy_arr, masked,
vmax=225., vmin=75., edgecolors='none'))
axarr[3].set_title('COSMIC Topside Scale Height')
axarr[3].set_ylabel('Solar Local Time (Hours)')
cbar = fig.colorbar(cax[3], ax=axarr[3])
cbar.set_label('Scale Height (km)')
axarr[3].set_xlabel('Apex Longitude (Degrees)')
fig.tight_layout()
fig.savefig('ssnl_median_ivm_cosmic_1d.png')
# Make COSMIC profile plots
# 6 pages of plots, 4 plots per page
for k in np.arange(6):
fig, axarr = plt.subplots(4, sharex=True, figsize=(8.5, 11))
# Iterate over a group of four sectors at a time (4 plots per page)
for (j, sector) in enumerate(list(
zip(*cosmic_results['ELEC_dens']['median']))[k * 4:(k + 1) * 4]):
# Iterate over all local times within longitude sector.
# Data is returned from the median routine in plot order, [y, x]
# instead of [x,y].
for (i, ltview) in enumerate(sector):
if np.isfinite(ltview):
# Plot a given longitude/local time profile
temp = pds.DataFrame(ltview['ELEC_dens'])
# Produce a grid covering plot region
# (y values determined by profile)
xx, yy = np.meshgrid(np.array([i, i + 1]),
np.arange(len(temp.index.values) + 1))
filtered = ma.array(np.log10(temp.values),
mask=pds.isnull(temp))
graph = axarr[j].pcolormesh(xx, yy, filtered,
vmin=3., vmax=6.5)
cbar = fig.colorbar(graph, ax=axarr[j])
cbar.set_label('Log Density (N/cc)')
axarr[j].set_xlim(0, 24)
axarr[j].set_ylim(0., 300.)
axarr[j].set_yticks([50., 100., 150., 200., 250.],
[150., 300., 450., 600., 750.])
axarr[j].set_ylabel('Altitude (km)')
axarr[j].set_title('Apex Longitudes %i-%i' %
(4 * k * 15 + j * 15, 4 * k * 15 + (j + 1) * 15))
axarr[-1].set_xticks([0., 6., 12., 18., 24.])
axarr[-1].set_xlabel('Solar Local Time of Profile Maximum Density')
fig.tight_layout()
fig.savefig('cosmic_part{}.png'.format(k))