.. _ex-general:
Load a general data file (DMSP SSIES)
=====================================
DMSP SSIES provides commonly used polar data, which can be accessed from
`Madrigal `_, which also has a Python API called
`madrigalWeb `_. To run this example,
follow the previous link(s) and download the ASCII file for F15 on 23 May 2000.
Choosing the UT DMSP with quality flags for the best calcuations of ion drift
and selecting ASCII instead of HDF5 will provide you with a file named
**dms_ut_20000523_15.002.txt**. To load this file, use the following commands.
::
import datetime as dt
import numpy as np
import matplotlib as mlt
import matplotlib.pyplot as plt
import ocbpy
hh = ["YEAR MONTH DAY HOUR MIN SEC RECNO KINDAT KINST UT1_UNIX UT2_UNIX GDALT GDLAT GLON MLAT MLT ION_V_SAT_FOR ION_V_SAT_LEFT VERT_ION_V NI PO+ PHE+ PH+ TI TE RPA_FLAG_UT IDM_FLAG_UT RMS_X SIGMA_VY SIGMA_VZ"]
dmsp_filename = "dms_ut_20000423_15.002.txt"
dmsp_head, dmsp_data = ocbpy.instruments.general.load_ascii_data(dmsp_filename, 1, datetime_cols=[0, 1, 2, 3, 4, 5], header=hh, datetime_fmt="%Y %m %d %H %M %S", int_cols=[6, 7, 8, 25, 26])
print(dmsp_data['TI'].shape, dmsp_data.keys())
(200060,), ['PH+', 'KINST', 'MIN', 'RPA_FLAG_UT', 'KINDAT', 'datetime', 'MLAT', 'UT2_UNIX', 'ION_V_SAT_FOR', 'ION_V_SAT_LEFT', 'GDALT', 'UT1_UNIX', 'GDLAT', 'HOUR', 'PHE+', 'IDM_FLAG_UT', 'SIGMA_VZ', 'SIGMA_VY', 'SEC', 'RMS_X', 'TI', 'TE', 'DAY', 'GLON', 'NI', 'RECNO', 'PO+', 'MLT', 'YEAR', 'MONTH', 'VERT_ION_V']
In the call to :py:func:`ocbpy.instruments.general.load_ascii_data`, quality
flags and number of points are saved as integers by specifying int_cols. The
header may not need to be specified using :py:data:`header`; have a go loading
it without this keyword arguement. The results will be the same as long as
there are no errors in the file header specification.
Before calculating the OCB coordinates, add space in the data dictionary for the
OCB coordinates and find out which data have a good quality flag.
::
ram_key = 'ION_V_SAT_FOR'
rpa_key = 'RPA_FLAG_UT'
idm_key = 'IDM_FLAG_UT'
dmsp_data['OCB_MLT'] = np.full(shape=dmsp_data[ram_key].shape,
fill_value=np.nan)
dmsp_data['OCB_LAT'] = np.full(shape=dmsp_data[ram_key].shape,
fill_value=np.nan)
igood = np.where((dmsp_data[rpa_key] < 3) & (dmsp_data[idm_key] < 3))
print(len(igood[0]), dmsp_data[ram_key][igood].max(),
dmsp_data[ram_key][igood].min())
7623 978.0 -2159.0
Now get the OCB coordinates for each location. This will not be possible
everywhere, since IMAGE doesn't provide Southern Hemisphere data and only times
with a good OCB established within the last 5 minutes will be used.
::
idmsp = 0
ndmsp = len(igood[0])
ocb = ocbpy.OCBoundary()
ocb.get_next_good_ocb_ind()
print(idmsp, ndmsp, ocb.rec_ind, ocb.records)
0 7623 0 305805
This is the starting point for cycling through the records.
::
while idmsp < ndmsp and ocb.rec_ind < ocb.records:
idmsp = ocbpy.match_data_ocb(ocb, dmsp_data['datetime'][igood],
idat=idmsp, max_tol=600)
if idmsp < ndmsp and ocb.rec_ind < ocb.records:
nlat, nmlt, r_corr = ocb.normal_coord(
dmsp_data['MLAT'][igood[0][idmsp]],
dmsp_data['MLT'][igood[0][idmsp]])
dmsp_data['OCB_LAT'][igood[0][idmsp]] = nlat
dmsp_data['OCB_MLT'][igood[0][idmsp]] = nmlt
idmsp += 1
igood = np.where(~np.isnan(dmsp_data['OCB_LAT']))
print(len(igood[0]), dmsp_data['OCB_LAT'][igood].max())
2218 89.5098551674689
Now, let's plot the satellite track over the pole, relative to the OCB, with
the location accouting for changes in the OCB at a 5 minute resolution. Note
how the resolution results in apparent jumps in the satellite location. We
aren't going to plot the ion velocity here, because it is provided in spacecraft
coordinates rather than magnetic coordinates, adding an additional
(and not intensive) level of processing.
::
fig = plt.figure()
fig.suptitle("DMSP F15 in OCB Coordinates on {:}".format(
dmsp_data['datetime'][igood][0].strftime('%d %B %Y')))
ax = fig.add_subplot(111, projection="polar")
ax.set_theta_zero_location("S")
ax.xaxis.set_ticks([0, 0.5 * np.pi, np.pi, 1.5 * np.pi])
ax.xaxis.set_ticklabels(["00:00", "06:00", "12:00 MLT", "18:00"])
ax.set_rlim(0, 40)
ax.set_rticks([10, 20, 30, 40])
ax.yaxis.set_ticklabels(["80$^\circ$", "70$^\circ$", "60$^\circ$",
"50$^\circ$"])
lon = np.arange(0.0, 2.0 * np.pi + 0.1, 0.1)
lat = np.ones(shape=lon.shape) * (90.0 - ocb.boundary_lat)
ax.plot(lon, lat, "m-", linewidth=2, label="OCB")
dmsp_lon = dmsp_data['OCB_MLT'][igood] * np.pi / 12.0
dmsp_lat = 90.0 - dmsp_data['OCB_LAT'][igood]
dmsp_time = mpl.dates.date2num(dmsp_data['datetime'][igood])
ax.scatter(dmsp_lon, dmsp_lat, c=dmsp_time, cmap=mpl.cm.get_cmap("viridis"),
marker="o", s=10)
ax.text(10 * np.pi / 12.0, 41, "Start of satellite track")
tticks = np.linspace(dmsp_time.min(), dmsp_time.max(), 6, endpoint=True)
dticks = ["{:s}".format(mpl.dates.num2date(tval).strftime("%H:%M"))
for tval in tticks]
cb = fig.colorbar(ax.collections[0], ax=ax, ticks=tticks,
orientation='horizontal')
cb.ax.set_xticklabels(dticks)
cb.set_label('Universal Time (HH:MM)')
ax.legend(fontsize='medium', bbox_to_anchor=(0.0,1.0))
.. image:: ../figures/example_dmsp_north_location.png