Classification of High Density Regions in Global Ionospheric Maps With Neural Networks

The database of Global Ionospheric Maps (GIMs) produced at Jet Propulsion Laboratory is analyzed. We define high density total electron content (TEC) regions (HDRs) in a map, following certain selection criteria. For the first time, we trained four convolutional neural networks (CNNs) corresponding to four phases of a solar cycle to classify the GIMs by the number of HDRs in each map with ∼ 80% accuracy on average. We compared HDR counts for GIMs across ten years to draw conclusions on how the number of HDRs in the GIMs changes throughout the solar cycle. Occurrence of HDRs during different geomagnetic activity conditions is discussed. Catalog of selected HDRs for ten years and four CNN‐based models that can be used to extend classification to other years are provided for the community to use.


Introduction
Ionosphere is a part of the Earth's upper atmosphere starting from the altitude of about 100 km up that is, partially ionized by solar EUV radiation. It responds promptly to disturbances in the Sun-Earth environment or space weather events (Eddy, 2009;Kelley, 1989). Abundance of electrons in the ionosphere is characterized by the total electron content (TEC), that is, the columnar number density of electrons (the total number of electrons obtained by integrating along a vertical tube with a cross section of one meter squared). The TEC is measured in total electron content units (TECU = 10 16 electrons/m 2 ). 2D distribution of TEC over the globe is visualized by Global Ionospheric Maps (GIMs, see the next section for the methodology), for example, (Mannucci et al., 1998;Roma-Dollase et al., 2018).
Ionospheric dynamics is controlled by solar and Earth's magnetospheric influences from "above," and by lower and middle atmospheric influences (tides, upward propagating atmospheric waves) from "below." Neutral counterpart of the upper atmosphere or thermosphere is closely coupled with the ionosphere through ion-neutral collisions and heating (Kelley, 1989;Schunk & Nagy, 2009). Equatorial Ionization Anomalies (EIAs) are the prominent feature of daytime ionosphere that are caused by the combined effect of direct solar influences and complex processes in the coupled magnetosphere-ionosphere-thermosphere (MIT) system or "equatorial fountain" (Appleton, 1954;Kelley, 1989). EIAs are two extended regions of elevated TEC located on opposite sides of geomagnetic equator (see Figure 1a). However, sometimes three distinct anomalies (Astafyeva et al., 2016(Astafyeva et al., , 2017Maruayama et al., 2016) or more are observed (see Figures 1b  and 1c). These multiple daytime TEC intensifications may not fit into strict definition of EIAs. However, formation of such dense ionization regions could provide important clues to dynamics of the MIT system and underlying physical processes in the ionosphere. We believe there is need for a systematic statistical study of EIA-like features in daytime ionosphere. However, a database of such ionospheric features, their statistics and analysis of connection to the solar cycle and geomagnetic activity are missing. The goal of this paper is to present a robust method to identify and classify high density TEC features and provide the first results of the data set analysis.
Artificial Intelligence (AI) and machine learning techniques have made possible new data analysis and discoveries in space physics and many other disciplines. Camporeale (2019) overviewed promising applications of machine learning to nowcasting and forecasting of space weather. The textbook (Camporeale et al., 2018) for the first time summarized a variety of machine learning methods and results for specific space weather problems. McGranaghan et al. (2017) demonstrated successful application of network analysis to an ionospheric data set and discovered connections between multi-scale TEC properties and solar influences. Machine learning shows promise in TEC forecasting (Liu et al., 2020). Chen et al. (2019) Abstract The database of Global Ionospheric Maps (GIMs) produced at Jet Propulsion Laboratory is analyzed. We define high density total electron content (TEC) regions (HDRs) in a map, following certain selection criteria. For the first time, we trained four convolutional neural networks (CNNs) corresponding to four phases of a solar cycle to classify the GIMs by the number of HDRs in each map with ∼80% accuracy on average. We compared HDR counts for GIMs across ten years to draw conclusions on how the number of HDRs in the GIMs changes throughout the solar cycle. Occurrence of HDRs during different geomagnetic activity conditions is discussed. Catalog of selected HDRs for ten years and four CNN-based models that can be used to extend classification to other years are provided for the community to use.
VERKHOGLYADOVA ET AL. demonstrated a deep learning based approach to fill in missing information in TEC maps. Convolution neural networks (CNNs) are deep learning AI networks which make use of a series of convolutional and pooling layers to extract important features from data (Goodfellow et al., 2016). CNNs have become popular for image classification problems because they are able to identify and extract important features in an image in order to differentiate between images of different objects. For example, a simple CNN could be trained to determine if any given image of a vehicle is a car, truck, or motorcycle. The CNN would first have to be "trained" using a few thousand labeled images of vehicles. During this "training," the CNN is given VERKHOGLYADOVA ET AL.  We identify two, three and five TEC intensifications, correspondingly, by eye. Ground sites used for GIM reconstruction each day are shown by magenta dots. example images along with the desired classification/output. As it is given more and more of these labeled "training" images, the network uses a process called "back propagation" to "learn" what combination of image features gives it the highest probability of correct vehicle classification. CNNs abilities for object recognition and image classification problems yet to be applied to the wealth of Earth's ionospheric data. We apply CNN approach to several years of GIM data and provide the first statistical analysis of solar cycle and geomagnetic activity dependencies of high density TEC features. The training data sets and models for automated classification are made publicly available as a result of this study.
We introduce the GIM data set in the next section. In the following sections we describe CNN-based approach used in the paper, discuss the method limitations and analyze GIM classification results. We will analyze occurrences of multiple HDRs per map in connection with a solar cycle phase and geomagnetic activity. Implications for understanding physical processes in daytime ionosphere will be discussed. Our main results are summarized in the Conclusions.

Data Set
GPS measurements provide integrated electron density along the ray path between a GPS satellite and a ground receiver. Slant TEC is typically inferred from delays of electromagnetic wave propagation through the ionosphere at L1 (1,575.42 MHz) and L2 (1,227.60 MHz) frequencies. There are several approaches to convert slant TEC to vertical TEC and create GIMs, see reviews by Hernandez-Pajares et al. (2017); Roma-Dollase et al. (2018). In this study we call vertical TEC as simply TEC. The TEC is measured in TECU = 10 16 electrons/m 2 ). We utilize the publicly available data set (https://sideshow.jpl.nasa.gov/pub/iono_daily/ gim_for_research/jpli/) based on TEC mapping technique developed at Jet Propulsion Laboratory (Iijima et al., 1999;Komjathy et al., 2005;Mannucci et al., 1998). The data set contains gridded TEC data on 1° by 1° horizontal grid generated using measurements from about 200 ground sites every 15 min, currently from 2000 to 2017. It should be noted that although GIMs (created at JPL or other institutions) are widely used to analyze global ionospheric dynamics, they have certain limitations. First, a limited number of stations are used for creating a map and JPL processing utilizes only ground-based sites with sparse distribution over oceans. Second, conversion of slant TEC to vertical TEC and subsequent interpolation over regions without TEC measurements adds uncertainty to GIM as a gridded TEC data product defined over the whole globe. Several evaluations of different GIM datasets, for example, (Roma-Dollase et al., 2018), using independent measurements showed that a typical error is below three TECU outside of high latitudes. TEC mapping technique is constantly evolving. Recently, Hernandez-Pajares et al. (2020) introduced a new approach to improve GIMs by utilizing measurements gathered from multi-platform spacecraft and oceanic vessels. An example of GIM map from the JPL database is shown in Figure 1. Magenta dots indicate ground sites used to create the maps. Different sets of ground sites can be used for creation of different maps, depending on data availability and quality. JPL algorithm selects an optimal distribution of available sites for a GIM.
We define high density TEC regions (HDRs) as regions on a map that satisfy the following criteria 1. TEC value of at least 25 TECU. 2. Significantly higher TEC than the surrounding region of the map. 3. Minimum of three GPS stations located within the region.
The first criterion limits our analysis to daytime TEC and at the same time excludes relatively faint TEC intensifications from being counted. The drawback of the limitation is that zero HDRs are identified for some of the maps during solar minimum years. We suggest that this is the first step in approach to select the most prominent TEC features and it is somewhat biased toward ionosphere near solar maximum years. Such threshold is helpful to define a "noise" level for a map. We will be experimenting with threshold selection based on individual characteristics of the map and taking into account solar cycle phase and geomagnetic activity. Following second criterion, we aim to discard faint local TEC intensifications. TEC gradient at HDR boundary will be calculated to ensure robust classification of a distinct feature (see the next section). The threshold for a distinct TEC intensification above a background is determined experimentally. The third criterion ensures that actual measurements within a high density region were utilized. Since GIM generation procedure is based on interpolation, this condition helps to reduce TEC mapping errors and increase reliability of HDR identification.
The solar cycle is one of the main factors that impacts TEC in the ionosphere. Generally, we can separate the cyclic solar cycle progression into four main phases: the solar minimum, the solar inclining phase, the solar maximum, and the solar declining phase. We will quantify the main features of the GIMs and explore their occurrence rates throughout the solar cycle.

Method
GIM data set of 18 years contains over 630,000 individual maps. It is not feasible to classify the maps by eye. Our approach to automated classification is two-fold. First, we perform "labeling" of a large set of maps to create training and testing datasets for use in CNN development. We utilize image processing routines described below to perform initial classification and verify the results by visual inspection. To make GIM classification more efficient, we train CNNs to perform classification of an extended data set. These steps are outlined below.
We create labeled data set by sorting GIMs into classes with different numbers of HDRs. Size of a region where TEC intensification is identified is not predetermined and an HDR is selected as long as three criteria are satisfied. In order to automate the data labeling process, we made use of OpenCV's image processing library for Python (Open Source Computer Vision (OpenCV), 2020). The Laplacian operator is often used to detect edges in an image. The operator is highly sensitive to noise in an image (Wang, 2007). We used the Laplacian operator from the OpenCV to calculate the Laplacian for each point on each map and set a minimum absolute-value threshold of −0.4 for the Laplacian value in order to focus on regions with significantly higher TEC than the surrounding area. Wide range of thresholds were tested and classification results were qualitatively compared to find an optimal value. It was found that when using threshold values smaller than −0.5, many areas which appear to be HDRs tended to be ignored. When threshold values larger than −0.3 were chosen, erroneous HDRs were identified. Next, we used the Dilate () and Erode () methods (Eroding and Dilating -OpenCV, 2019) to group together these dense regions into individual objects which could be identified by openCV's Connected Components operator. This operator identifies and counts all individual objects in an image (OpenCV Structural Analysis and Shape Descriptors, 2020). Figure 2 illustrates the image processing steps to classify HDR in a GIM image.
We conducted experiments to test the robustness of using the Laplacian operator to identify HDRs in images. We found out that labeling script works successfully when there are visually identified sharp edges to a TEC brightening. Otherwise, the Laplacian map can be noisy and would impede correct count of HDRs in a map. GIMs are the result of interpolation of measurements from an optimally selected set of ground VERKHOGLYADOVA ET AL.  sites. These sites are unevenly distributed over the globe with preferential locations over continents. Since interpolation commonly results in smoothing over sharp gradients rather than creating unphysical ones, we consider our classification method as conservative. Since our noise floor is 25 TECU, we avoid selecting potential unphysical brightenings over areas with sparse site distributions. We are fairly confident in selection of an HDR as relatively bright TEC region (by the TEC magnitude) compared to neighboring regions. Since success of the automated labeling relies partially on sharpness of the main features of a TEC map, visual inspection is necessary to correctly identify HDRs. This consideration prompted us to use trained CNNs to classify a large data set of GIMs rather than rely on labeling script alone.
Using this selection process, we were able to generate estimations for the numbers of HDRs in each GIM. Maps where the Connected Components function identified regions which were very close together (indicating that the two HDRs identified might actually be one large HDR), unusually large (indicating that the large HDR identified might actually be several smaller HDRs), or were too small or misshapen, were visually inspected to confirm the number of HDRs. We found that the labeling script is helpful to aid selection of HDRs but alone is not sufficient for successful map classification. On the other hand, visual inspection and defining of the number of HDRs on a map by eye is a subjective approach and could introduce individual biases in classification of complex TEC intensifications. Examples are identifications of three and five HDRs in Figure 1b and Figure 1c. Smooth transitions between different HDRs in some GIMs present limitations for our image-based classification technique. We are exploring using additional methods to sharpen images. In our best effort current approach, we utilize labeling algorithm in combination with three criteria and visual inspection to classify GIMs for training and testing CNNs.
CNN approach has the means to automate image classification process and determine the morphological types of larger quantities of images (Goodfellow et al., 2016). To classify GIMs, we utilize an open-source Python API Keras (https://keras.io). It allows for easy front-end design of many different types of deep AI networks, including CNNs. Keras makes available many pre-trained CNNs which have been trained on large data sets to extract basic image features for image classification. Users can utilize the pre-trained CNNs and train them using a set of images for their specific image classification problem. For example, one of these networks could be trained using a relatively small set of labeled images in order to create a CNN to look for and sort by a certain feature, for example, isolated bright regions. The benefit of using such pretrained networks is that, since they are already equipped to extract basic image features, a smaller training data set can be used and a high degree of accuracy for a specific classification problem can be still obtained. One such pre-trained network is MobileNet available from Keras which has become increasingly popular for image classification (Howard et al., 2017).
Our second step was to train a MobileNet. We labeled approximately 10,000 TEC maps, taken from ten different years in the JPL database (2000, 2001, 2003, 2004, 2008, 2009, 2010, 2011, 2014, and 2015), with the number of HDRs. Of these ∼10,000 labeled maps, ∼9,000 were used to train MobileNet to classify the maps by the number of HDRs. The remaining ∼1,000 labeled images, which we will refer to as the "test set," were set aside to test the accuracy of the CNN after training. Approximately an equal number of training maps were taken from each year and same for testing maps. All seasons were sampled and no preference was given to UT times of GIMs. Two and three HDRs appear more frequently in the labeled data set, followed by occurrences of four HDRs.
After training MobileNet using the "training set" images, we found that the resultant trained CNN was approximately 76% successful in correctly classifying the test set images, which is not optimal if we want to be able to rely on this network to predict the number of HDRs in other TEC maps accurately. Table 1 lists the percentage of successfully identified HDRs per the HDR number. The lower numbers of HDRs per map were the most successfully classified. It is possible that the CNN did a better job of predicting test set images with numbers of HDRs that were more frequency represented in the training set (two and three HDRs). We VERKHOGLYADOVA ET AL. also expect that the case of two HDRs should be the most typical. Note an increase in successful identification of five through six HDRs.
The majority of cases where the network predicted the incorrect number of HDRs were cases where the network was only off by one HDR count rather than two or three. We found that the network can reliably predict the number of HDRs to within one of the actual number over 98% of the time. Using the trained CNN, we were able to fully automate classification of the number of HDRs in the remaining unlabeled TEC maps in each of the 10 years (2000, 2001, 2003, 2004, 2008, 2009, 2010, 2011, 2014, and 2015) from the JPL GIM database.

Solar EUV radiation controls ionization of the upper Earth's atmosphere and affects daytime TEC values.
Since the EUV radiation varies throughout a solar cycle, the resulting global TEC patterns and values depend on a solar cycle phase. We tested the hypothesis that we get better results by training four specialized networks for each particular solar cycle phase: solar inclining phase, solar maximum, solar declining phase and solar minimum. In order to improve the accuracy of our network, we divided up the training data into four groups representing the four different phases in a solar cycle.  Figure 3. By separating our data into these four groups, we then trained four separate networks, one network for each solar cycle phase. Figure 4 provides a schematic diagram of the process to select and label data, train and test CNNs.
The four individual specialized networks we trained were more accurate at ∼80% in classifying the test set of images, but they were all still below 85% accuracy. For example, CNN for solar minimum demonstarted 82% accuracy. Using a larger training data set could improve the accuracy of our network (Jain & Chandrasekaran, 1982;Raudys & Jain, 1991). We noticed that each of these networks is able to predict within one HDR count of the correct number over 95% of the time. Detailed analysis of HDR classification sensitivity to the threshold value (selected here as −0.4) warrants a separate study. These trained networks are made available online for further experimenting and testing.

Analysis of Classification Results for Global Ionospheric Maps
We used these four specialized CNNs to determine a prediction for the number of HDRs in all of the TEC maps for ten years of data from the JPL database. For example, since 2008 is a solar minimum year, we utilized the solar minimum CNN to classify all of the TEC maps from 2008 and determine the number of VERKHOGLYADOVA ET AL.
10.1029/2021EA001639 6 of 12 HDRs for each map. We then used these classifications to calculate the relative frequency of TEC maps with 0, 1, 2, 3, 4, 5, and 6 HDRs for each year. We found that only approximately 37% of the TEC maps from 2008 had two HDRs. The remaining two-thirds of classified GIM had either very faint TEC intensifications (classified as zero or one HDR) or multiple high-density TEC regions were identified.
In order to more clearly illustrate the relationship between HDR frequencies and the solar cycle, we plotted aggregate statistics of the percentage of HDRs with different counts for years in the same solar cycle phase (see Figure 5). For example, the blue solar minimum bars show the average of the frequencies we calculated for 2008 and 2009 which are the two solar minimum years.
Looking at this Figure 5, we see the following trends. The solar minimum years (shown in blue) generally tend to have the fewest numbers of HDRs, peaking at two HDRs per map, while the solar maximum years, VERKHOGLYADOVA ET AL.
10.1029/2021EA001639 7 of 12  which are shown in gray, tend to have larger numbers of HDRs, peaking around three or four HDRs per map. The solar inclining and declining years seem to fall in between. The large number of zero or single HDRs found in solar minimum years indicate a potential shortcoming of our approach. Lowering the TEC threshold from 25 TECU could result in increase of the number of identified HDRs during low solar activity times. However, lowering acceptable noise level would also result in identification of smaller TEC intensifications that may well be within TEC uncertainty of mapping (especially in areas with sparse distribution of ground sites). We are currently researching optimal threshold identification for individual maps.
The multiple HDR occurrences in solar maximum years (and to a lesser extent during solar declining phase) are consistent with the view that increased solar activity causes larger daytime TEC. That could result in more high density TEC regions. We would like to note that zero (corresponding to low values of TEC in a map), one and two HDRs occur more frequently during solar inclining phase years as compared to the solar declining phase years. Three to five HRDS occur more frequently during solar declining phase years as compared to the solar inclining phase years. One possible explanation is that more intense geomagnetic storms tend to occur during declining phases than inclining phases, as indicated by several studies (Bell et al., 1997;Gonzalez et al., 2011), if more multiple-HDRs appear during storm time than quiet time. We will explore a dependence of HDR counts on geomagnetic activity below. Furthermore, it is possible that the number of HDRs are not only driven by the solar cycle directly but also by interaction with middle and low atmosphere through tides and large-scale waves. Thus, the solar EUV cycle could be the main, but not the only factor defining occurrence of multiple HDRs. Astafyeva et al. (2016) and Maruayama et al. (2016) analyzed occurrence of TEC anomalies in addition to EIAs or three-peak EIA structures. They attributed modification of the fountain mechanism to thermospheric winds. Enhanced neutral winds during geomagnetically quiet times drive generation of threepeak structures predominantly during hemispheric summers (May-August in the northern hemisphere and November-January in the Southern hemisphere) and close to solstices. Detailed statistical analysis for 2003-2009 years revealed increased occurrences during 2005, 2008, and 2009 years, that correspond to solar minimum and declining phases (Astafyeva et al., 2016). Our results shown in Figure 5 show a different trend with the largest counts of three HDRs per map occurring during solar maximum and declining phase years. However, HDRs in this study are defined more broadly than EIAs (or three-peak structures of EIAs) as any distinct high-density regions and include EIAs as a subset.
Further detailed analysis is performed on the training (or labeled) data set that samples TEC over ten years. We estimate the cumulative number of different global HDR counts (1, 2, 3, and 4) for each day of a year (DOY). This allows us to trace seasonal preferences if any. Figure 6 shows moving averages with 27-days window of the total number of different HDR counts for each DOY. We chose the running averages window to be equal to the synodic solar rotation period that corresponds to a rotation period of a fixed feature, for example, active region, on the Sun. This smoothing accounts for the known modulation of the low-latitude ionospheric TEC (Min et al., 2009) and reduces its effect on our results. A single HDR over the globe occurs with the least frequency over the year as compared to multiple HDR occurrences. It corresponds to relatively weak EIAs that are difficult to distinguish into two separate TEC intensification regions. They appear more frequently during summer and winter months. Our analysis does not allow to distinguish between appearances of HDRs in different hemispheres. Such additional analysis is particularly difficult if an HDR spans two hemispheres, unlike cases of three-peak EIAs where the third peak is clearly visible in one or another hemisphere. The counts of two, three and four distinct HDRs appear more frequently during spring and autumn seasons. A local maximum is notable for three-peak structures around 150-201 DOY (approximately June through August) that may correspond to the statistical results of (Astafyeva et al., 2016) for three-peak EIAs. Three HDRs occur the most frequently around 70-90 DOY. Several local maxima are observed for occurrences of three and four HDRs. Two and three HDRs in a map occur more frequently than one or four HDRs. Note that the cases of three HDRs may include EIAs with complex structure inside that our algorithm identified as separate density regions and three-peak structures. The occurrences of three HDRs are about equal to occurrences of two HDRs that include typical two-peak EIAs.
To address dependence of HDR counts on geomagnetic activity, we analyze distribution of the counts with 3 h Kp index. Figure 7a shows that the dominant configuration of a GIM at low values of Kp index is a two-peak TEC enhancement that likely corresponds to EIAs. Three-peak HDRs occur next frequently and predominantly at low levels of geomagnetic activity. Three and more HDRs per map are also observed at moderately disturbed conditions with 3 < Kp < 4. Note, that were are very few instances (less than 100 occurrences) over the 10-years training set selection with large Kp > 4 values. Hourly Dst index indicates geomagnetic storm activity in low-latitude regions. The histogram (see Figure 7b) shows dominance of twopeak and three-peak structures during weak geomagnetic activity (Dst > −20 nT). The data set includes storm intervals with Dst around −400 nT but they contribute a small number of HDR counts compared with moderate storm intervals (Dst > −100 nT). Single HDR in a GIM map occurs during quiet periods.

Conclusion
For the first time, we applied transfer learning with MobileNet to classify the number of HDRs in global distribution of TEC over ten years (solar cycles 23 and 24). The JPL database of GIMs was utilized in the study. We applied simple HDR selection criteria and found that using four individual CNNs trained for the phases of a solar cycle achieved the HDR classification accuracy of ∼80% on average.
There is a pronounced dependence of how frequently a certain number of HDRs occurs in a map on a solar cycle phase. We found that two and three HDRs are the dominant features of GIMs throughout different phases of a solar cycle. Two HDRs per map likely correspond to EIAs. However, occurrences of multiple HDRs in a map indicate that in many cases EIAs either have a complex inhomogeneous structure or secondary high-density TEC regions are present. Figure 5 indicates that three HDRs per map appear slightly more frequently than two HDRs. Larger numbers of HDRs per map (from three to five) occur more frequently during solar maximum and declining phase. GIMs during solar minimum years tend to have two HDRs or less. Close to a solar minimum in 2008, only ∼37% of the GIMs in have two HDRs, while the most of the maps have either zero or one HDRs. A single HDR or no HDR per map indicates that TEC enhancements are too weak to be identified by our approach.
The counts of two, three and four distinct HDRs appear more frequently during spring and autumn seasons. Three-peak EIAs also tend to occur around June through August and are potentially a subset of three HDRs occurrences. This result is in a general agreement with prior analysis of three-peak EIAs by Maruayama et al. (2016) and Astafyeva et al. (2016). It is possible that thermospheric wind anomalies are responsible for modification of the fountain mechanism leading to creation of multiple HDRs. Our analysis does not show VERKHOGLYADOVA ET AL.
10.1029/2021EA001639 9 of 12 a preference for multiple HDRs to occur during periods of elevated geomagnetic activity. However, there is a secondary peak in occurrences of three HDRs per map at Kp > 3. Understanding frequent occurrences of multiple dense TEC regions in global maps warrants further study. Elucidation of a physical mechanism responsible for formation of multiple HDRs requires further observational analysis and modeling efforts, see for example, (Astafyeva et al., 2017).
For the first time we created four CNN-based models that can be utilized to classify GIMs and catalog HDRs for future analysis by the community. We realize limitations of the methodology but these initial VERKHOGLYADOVA ET AL.
10.1029/2021EA001639 10 of 12 classification results can help to analyze and explore under what conditions multiple large-scale TEC intensification can occur globally. This multi-year global data set can aid our understanding of large-scale structuring of the daytime ionosphere under different driving conditions.

Data Availability Statement
The catalog of HDR counts for GIMs in the training and testing sets, and four CNN models trained for solar cycle phases are provided in the online data set (https://doi.org/10.6084/m9.figshare.13501872) associated with the paper. The paper utilized GIM data set that is, provided at https://sideshow.jpl.nasa.gov/pub/ iono_daily/gim_for_research/jpli/.