Merge pull request #1342 from spedas/t_1341

Adds support to load USGS Variometer data, adds test cases for loading variometer data.

Author: jameswilburlewis

pyspedas/projects/themis/ground/gmag.py

@@ -44,6 +44,7 @@ def gmag(
     no_update=False,
     time_clip=False,
     force_download=False,
+    sampling_rate=1,
 ):
     """
     Load ground magnetometer data from the THEMIS mission.
@@ -84,7 +85,9 @@ def gmag(
     force_download: bool, optional
         Download file even if local version is more recent than server version
         Default: False
-
+    sampling_rate: int, optional
+        Specify a sampling rate for loading variometer data. Accepts 1 (Hz) or 10 (Hz).
+        Default: 1
     Returns
     -------
     dict
@@ -94,10 +97,28 @@ def gmag(
     --------
     >>> from pyspedas.projects.themis import gmag
     >>> from pyspedas import tplot
+    >>> from pyspedas import subtract_median
     >>>
     >>> # Load ground magnetometer data for specific sites and time range
     >>> gmag_vars = gmag(sites=['ccnv','bmls'], trange=['2013-11-05', '2013-11-06'])
     >>> tplot(['thg_mag_bmls', 'thg_mag_ccnv'])
+    >>>
+    >>> # Load variometer data for specific sites and time range:
+    >>> gmag_vars = gmag(sites=['s61a','anmo'], trange=['2026-02-24', '2026-02-25'])
+    >>> subtract_median(['thg_mag_s61a', 'thg_mag_anmo'])    
+    >>> tplot(['thg_mag_s61a-m', 'thg_mag_anmo-m'])
+    >>>
+    >>> # Load 10 Hz variometer data for specific sites and time range:
+    >>> gmag_vars = gmag(sites=['s61a','anmo'], sampling_rate=10, trange=['2026-02-24', '2026-02-25'])
+    >>> subtract_median(['thg_mag_s61a_100ms', 'thg_mag_anmo_100ms']) 
+    >>> tplot(['thg_mag_s61a_100ms-m', 'thg_mag_anmo_100ms-m'])
+    >>>
+    >>> # Load 10 Hz variometer data for specific sites and time range using the 10 Hz file name format:
+    >>> gmag_vars = gmag(sites=['s61a_100ms','anmo_100ms'], trange=['2026-02-24', '2026-02-25'])
+    >>> subtract_median(['thg_mag_s61a_100ms', 'thg_mag_anmo_100ms'])
+    >>> tplot(['thg_mag_s61a_100ms-m', 'thg_mag_anmo_100ms-m'])
+    >>>
+    
     """
 
     if sites is None:
@@ -250,6 +271,48 @@ def gmag(
             "weyb",
             "wgry",
         ]
+        variometer_sites=[
+            "anmo",
+            "casy",
+            "ccm",
+            "cola",
+            "cor",
+            "dgmt",
+            "dwpf",
+            "ecsd",
+            "eymn",
+            "e46a",
+            "e62a",
+            "goga",
+            "hrv",
+            "j47a",
+            "kbs",
+            "kevo",
+            "kono",
+            "ksu1",
+            "k30b",
+            "k50a",
+            "mbwa",
+            "mstx",
+            "m63a",
+            "o20a",
+            "pab",
+            "p57a",
+            "qspa",
+            "rssd",
+            "r49a",
+            "sba",
+            "sfjd",
+            "spmn",
+            "sspa",
+            "s61a",
+            "t47a",
+            "u38b",
+            "wci",
+            "whtx",
+            "wvt",
+            "352a"
+        ]
         sites = (
             thm_sites
             + tgo_sites
@@ -266,6 +329,7 @@ def gmag(
             + fmi_sites
             + aair_sites
             + carisma_sites
+            + variometer_sites
         )
 
     if group is not None:
@@ -274,6 +338,29 @@ def gmag(
     if not isinstance(sites, list):
         sites = [sites]
 
+    # check for sites in Variometer group
+    variometer = []
+    for site in sites:
+        s_rate=check_variometer(site)
+        # if check_variometer determines site
+        # to be 1 Hz sampling rate, check if 
+        # sampling_rate argument has been set
+        if (s_rate == 1) and (sampling_rate is not None):
+            # if so, pass sampling_rate value
+            # (case where variometer base names
+            # are passed, but 10 Hz sampling rate
+            # data is to be loaded)
+            
+            # check if sampling_rate is valid:
+            if sampling_rate not in [1,10]:
+                # throw error
+                logging.error('sampling_rate must be either 1 or 10. Setting to default value of 1')
+                sampling_rate=1
+            
+            variometer.append(sampling_rate) 
+        else:
+            variometer.append(s_rate)
+
     # check for sites in Greenland
     greenland = []
     for site in sites:
@@ -295,6 +382,7 @@ def gmag(
         notplot=notplot,
         stations=sites,
         greenland=greenland,
+        variometer=variometer,
         time_clip=time_clip,
         no_update=no_update,
         force_download=force_download,
@@ -482,3 +570,51 @@ def check_greenland(station_name):
                 return 1
 
     return 0
+
+def check_variometer(station_name):
+    """
+    Check if a given station belongs to the Variometers group. 
+    If so, check if the given stationname is formatted to request 
+    data in the 10 Hz sampling rate.
+
+    Parameters
+    ----------
+    station_name : str
+        The name of the station to check.
+
+    Returns
+    -------
+    int
+        Returns sampling rate (1 or 10) if the station is in the variometers group, otherwise returns 0.
+
+    Examples
+    --------
+    >>> from pyspedas.projects.themis.ground.gmag import check_variometer
+    >>> check_variometer("s61a")
+    1
+    >>> check_variometer("s61a_100ms")
+    10
+    """
+
+    # check if station_name contains "_100ms". If it does,
+    # create copy of station name with the "_100ms" 
+    # substring removed
+    if "_100ms" in station_name.lower():
+        station_name_base=station_name.lower()
+        station_name_base=station_name_base.replace("_100ms","")
+    else:
+        station_name_base=station_name.lower()
+    
+    gmag_dict = themis_gmag_dict.get_gmag_list()
+    for station in gmag_dict:
+        # compare a copy of station name with any "_100ms" 
+        # substring removed with the station ccode
+        if station["ccode"].lower() == station_name_base:
+            if (
+                station["variom"] is not None and station["variom"].lower() == "y"
+            ):
+                if "_100ms" in station_name.lower():
+                    return 10
+                else:
+                    return 1
+    return 0

pyspedas/projects/themis/load.py

@@ -14,6 +14,7 @@ def load(trange=['2013-11-5', '2013-11-6'],
          datatype=None, # ASK data, ESD (3d L2 ESA)
          stations=None,  # ground mag and ASK data
          greenland=None,  # also for ground mag data
+         variometer=None, # variometer check and sampling rate
          prefix='',
          suffix='',
          get_support_data=False,
@@ -177,10 +178,20 @@ def load(trange=['2013-11-5', '2013-11-6'],
                 return
             else:
                 pathformat = []
-                for site, in_greenland in zip(stations, greenland):
+                for site, in_greenland, s_variometer in zip(stations, greenland, variometer):
                     if site == 'idx':
                         # THEMIS GMAG index files are only L1
                         pathformat.append('thg/l1/mag/idx/%Y/thg_l1_idx_%Y%m%d_v??.cdf')
+                    elif s_variometer == 1:
+                        pathformat.append('thg/' + level + '/variometers/' + site
+                                          + '/%Y/thg_' + level + '_mag_' + site
+                                          + '_%Y%m%d_v??.cdf')
+                    elif s_variometer == 10:
+                        if "_100ms" in site.lower():
+                            site=site.replace("_100ms","")
+                        pathformat.append('thg/' + level + '/variometers/' + site
+                                          + '/%Y/thg_' + level + '_mag_' + site
+                                          + '_100ms_%Y%m%d_v??.cdf')
                     elif in_greenland:
                         pathformat.append('thg/greenland_gmag/' + level
                                           + '/' + site + '/%Y/thg_' + level

pyspedas/projects/themis/tests/test_themis.py

@@ -37,11 +37,40 @@ def test_check_gmag(self):
         self.assertTrue(check_gmag('ccnv') == 1)
         self.assertTrue(check_gmag('abcd') == 0)
 
+    def test_check_greenland(self):
+        """Check if a gmag station is in the greenland group."""
+        from pyspedas.projects.themis.ground.gmag import check_greenland
+        self.assertTrue(check_greenland('bfe') == 1)
+        self.assertTrue(check_greenland('bou') == 0)
+
+    def test_check_variometer(self):
+        """Check if a gmag station is in the 1 Hz variometer group or 10 Hz variometer group or not."""
+        from pyspedas.projects.themis.ground.gmag import check_variometer
+        self.assertTrue(check_variometer('anmo') == 1)
+        self.assertTrue(check_variometer('anmo_100ms') == 10)
+        self.assertTrue(check_variometer('bou') == 0)
+
     def test_load_gmag_data(self):
         """Load gmag."""
         pyspedas.projects.themis.gmag(varnames=['thg_mag_amer'], sites='amer')
         self.assertTrue(data_exists('thg_mag_amer'))
 
+    def test_load_gmag_variometer_1_hz_data(self):
+        """Load gmag variometer 1 Hz data."""
+        pyspedas.projects.themis.gmag(varnames=['thg_mag_s61a'], sites='s61a',trange=['2026-02-24', '2026-02-25'])
+        self.assertTrue(data_exists('thg_mag_s61a'))
+    
+    def test_load_gmag_variometer_10_hz_data_notag(self):
+        """Load gmag variometer 10 Hz data using the sampling_rate argument and not the time resolution tag."""
+        pyspedas.projects.themis.gmag(varnames=['thg_mag_s61a_100ms'], sites='s61a',sampling_rate=10,trange=['2026-02-24', '2026-02-25'])
+        self.assertTrue(data_exists('thg_mag_s61a_100ms'))
+
+    def test_load_gmag_variometer_10_hz_data_tag(self):
+        """Load gmag variometer 10 Hz data using the time resolution tag and not the sampling_rate argument."""
+        pyspedas.projects.themis.gmag(varnames=['thg_mag_s61a_100ms'], sites='s61a_100ms',trange=['2026-02-24', '2026-02-25'])
+        self.assertTrue(data_exists('thg_mag_s61a_100ms'))
+    
+
 
 class LoadTestCases(unittest.TestCase):
     """Test themis load functions."""