From dc9cc949873e649788e69a21dd835b8888b4fdce Mon Sep 17 00:00:00 2001
From: Nicholas Shavers <nshavers@contractor.usgs.gov>
Date: Thu, 5 Dec 2024 11:25:57 -0800
Subject: [PATCH] minor adjustments for spec. temperature channel
 compatibility.

---
 geomagio/Controller.py     |   5 ++
 geomagio/ImagCDFFactory.py | 112 +++++++++++++++----------------------
 geomagio/api/ws/Element.py |   3 +
 3 files changed, 53 insertions(+), 67 deletions(-)

diff --git a/geomagio/Controller.py b/geomagio/Controller.py
index 564b54c9..40e6b173 100644
--- a/geomagio/Controller.py
+++ b/geomagio/Controller.py
@@ -7,6 +7,7 @@ from typing import List, Optional, Tuple, Union
 
 from obspy.core import Stream, UTCDateTime
 
+from geomagio.ImagCDFFactory import ImagCDFFactory
 
 from .algorithm import Algorithm, algorithms, AlgorithmException, FilterAlgorithm
 from .DerivedTimeseriesFactory import DerivedTimeseriesFactory
@@ -27,6 +28,7 @@ from . import vbf
 from . import xml
 from . import covjson
 from . import netcdf
+from geomagio.ImagCDFFactory import ImagCDFFactory
 
 
 class Controller(object):
@@ -630,6 +632,8 @@ def get_output_factory(args):
         )
     elif output_type == "plot":
         output_factory = PlotTimeseriesFactory()
+    elif output_type == "imagcdf":
+        output_factory = ImagCDFFactory()
     else:
         # stream compatible factories
         if output_type == "binlog":
@@ -1040,6 +1044,7 @@ def parse_args(args):
             "xml",
             "covjson",
             "netcdf",
+            "imagcdf",
         ),
         # TODO: set default to 'iaga2002'
         help="Output format",
diff --git a/geomagio/ImagCDFFactory.py b/geomagio/ImagCDFFactory.py
index 3189d15a..9644c23e 100644
--- a/geomagio/ImagCDFFactory.py
+++ b/geomagio/ImagCDFFactory.py
@@ -6,8 +6,9 @@ The ImagCDF format is based on NASA's Common Data Format (CDF) and is designed
 to store geomagnetic data with high precision.
 
 References:
-- ImagCDF Format Documentation: https://intermagnet.org/docs/technical/im_tn_8_ImagCDF.pdf
-- CDF Library: http://cdf.gsfc.nasa.gov/
+- ImagCDF Technical Note: https://intermagnet.org/docs/technical/im_tn_8_ImagCDF.pdf
+- ImagCDF Official Documentation: https://tech-man.intermagnet.org/latest/appendices/dataformats.html#imagcdfv1-2-intermagnet-exchange-format
+- CDF: http://cdf.gsfc.nasa.gov/
 """
 
 from __future__ import absolute_import, print_function
@@ -20,6 +21,7 @@ import numpy as np
 from obspy import Stream, Trace, UTCDateTime
 
 from geomagio.TimeseriesFactory import TimeseriesFactory
+from geomagio.api.ws.Element import TEMPERATURE_ELEMENTS_ID
 
 from .geomag_types import DataInterval, DataType
 from .TimeseriesFactoryException import TimeseriesFactoryException
@@ -121,6 +123,7 @@ class ImagCDFFactory(TimeseriesFactory):
         Note: Parsing from strings is not implemented in this factory.
         """
         raise NotImplementedError('"parse_string" not implemented')
+    
     def write_file(self, fh, timeseries: Stream, channels: List[str]):
         # Create a temporary file to write the CDF data
         with tempfile.NamedTemporaryFile(delete=False) as tmp_file:
@@ -154,18 +157,24 @@ class ImagCDFFactory(TimeseriesFactory):
 
                 # Write time variable
                 cdf_writer.write_var(var_spec, var_attrs, ts_data)
-
+            
+            
             # Data variables
+            temperature_index = 0
             for trace in timeseries:
                 channel = trace.stats.channel
-                var_name = f"GeomagneticField{channel}"
+                if channel in TEMPERATURE_ELEMENTS_ID:
+                    temperature_index += 1 #MUST INCREMENT INDEX BEFORE USING
+                    var_name = f"Temperature{temperature_index}"
+                else:
+                    var_name = f"GeomagneticField{channel}"
                 data_type = self._get_cdf_data_type(trace)
                 num_elements = 1
                 CDF_CHAR = 51
                 CDF_UCHAR = 52
                 if data_type in [CDF_CHAR, CDF_UCHAR]:  # Handle string types
                     num_elements = len(trace.data[0]) if len(trace.data) > 0 else 1
-
+            
                 var_spec = {
                     'Variable': var_name,
                     'Data_Type': data_type,
@@ -178,7 +187,7 @@ class ImagCDFFactory(TimeseriesFactory):
                     'Pad': None,
                 }
 
-                var_attrs = self._create_var_attrs(trace)
+                var_attrs = self._create_var_attrs(trace, temperature_index)
 
                 # Write data variable
                 cdf_writer.write_var(var_spec, var_attrs, trace.data)
@@ -362,7 +371,8 @@ class ImagCDFFactory(TimeseriesFactory):
     def _create_time_stamp_variables(self, timeseries: Stream) -> dict:
         vector_times = None
         scalar_times = None
-
+        temperature_times = {}
+        temperature_index = 1
         for trace in timeseries:
             channel = trace.stats.channel
             times = [
@@ -371,7 +381,6 @@ class ImagCDFFactory(TimeseriesFactory):
             ]
             # Convert timestamps to TT2000 format required by CDF
             tt2000_times = cdflib.cdfepoch.timestamp_to_tt2000([time.timestamp() for time in times])
-            # tt2000_times = cdflib.cdfepoch.compute_tt2000(times) #this does not work
 
             if channel in self._get_vector_elements():
                 if vector_times is None:
@@ -385,16 +394,21 @@ class ImagCDFFactory(TimeseriesFactory):
                 else:
                     if not np.array_equal(scalar_times, tt2000_times):
                         raise ValueError("Time stamps for scalar channels are not the same.")
-            else:
-                # Handle other channels if necessary
-                pass
+            elif channel in TEMPERATURE_ELEMENTS_ID:
+                ts_key = f"Temperature{temperature_index}Times"
+                if ts_key not in temperature_times:
+                    temperature_times[ts_key] = tt2000_times
+                    temperature_index += 1
+                else:
+                    temperature_times[ts_key] = tt2000_times
 
         time_vars = {}
         if vector_times is not None:
             time_vars['GeomagneticVectorTimes'] = vector_times
         if scalar_times is not None:
             time_vars['GeomagneticScalarTimes'] = scalar_times
-
+        if temperature_times:
+            time_vars.update(temperature_times)
         return time_vars
 
 
@@ -441,8 +455,7 @@ class ImagCDFFactory(TimeseriesFactory):
         }
         return var_spec
 
-    def _create_var_attrs(self, trace: Trace) -> dict:
-        # print(trace.stats)
+    def _create_var_attrs(self, trace: Trace, temperature_index: Optional[int] = None) -> dict:
         channel = trace.stats.channel
         fieldnam = f"Geomagnetic Field Element {channel}" # “Geomagnetic Field Element ” + the element code or “Temperature ” + the name of the location where the temperature was recorded.
         units = '' # Must be one of “nT”, “Degrees of arc” or “Celsius”
@@ -454,10 +467,12 @@ class ImagCDFFactory(TimeseriesFactory):
             units = 'Degrees of arc'
             validmin = -90.0 
             validmax = 90.0 #The magnetic field vector can point straight down (+90°), horizontal (0°), or straight up (-90°).
-        elif 'Temperature' in channel:
+        elif channel in TEMPERATURE_ELEMENTS_ID:
             units = 'Celsius'
-            fieldnam = f"Temperature {trace.stats.location}"
-        elif channel == 'F':
+            fieldnam = f"Temperature {temperature_index} {trace.stats.location}"
+            validmin = -273.15
+            validmax = 79_999
+        elif channel in ['F','S']:
             units = 'nT'
             validmin = 0.0 # negative magnetic field intestity not physically meaningful.
             validmax = 79_999.0
@@ -466,12 +481,13 @@ class ImagCDFFactory(TimeseriesFactory):
             validmin = -79_999.0
             validmax = 79_999.0
 
+        # Determine DEPEND_0 based on channel type
         if channel in self._get_vector_elements():
             depend_0 = 'GeomagneticVectorTimes'
         elif channel in self._get_scalar_elements():
             depend_0 = 'GeomagneticScalarTimes'
-        else:
-            depend_0 = None  # Handle other cases if necessary
+        elif channel in TEMPERATURE_ELEMENTS_ID:
+            depend_0 = f"Temperature{temperature_index}Times"
 
     
         var_attrs = {
@@ -486,7 +502,6 @@ class ImagCDFFactory(TimeseriesFactory):
         }
         return var_attrs
 
-
     def _create_time_var_attrs(self, ts_name: str) -> dict:
         """
         Create a dictionary of time variable attributes.
@@ -546,7 +561,7 @@ class ImagCDFFactory(TimeseriesFactory):
         stream = Stream()
         # Read time variables
         time_vars = {}
-        for var in cdf.cdf_info()['zVariables']:
+        for var in cdf.cdf_info().zVariables:
             if var.endswith('Time'):
                 time_data = cdf.varget(var)
                 # Convert TT2000 to UTCDateTime
@@ -554,7 +569,7 @@ class ImagCDFFactory(TimeseriesFactory):
                 time_vars[var] = utc_times
 
         # Read data variables
-        for var in cdf.cdf_info()['zVariables']:
+        for var in cdf.cdf_info().zVariables:
             if not var.endswith('Time'):
                 data = cdf.varget(var)
                 attrs = cdf.varattsget(var)
@@ -635,14 +650,16 @@ class ImagCDFFactory(TimeseriesFactory):
 
         This method constructs the filename based on the ImagCDF naming
         conventions, which include the observatory code, date-time formatted
-        according to the data interval, and the publication level.
+        according to the data interval, and the publication level.  
+        
+        [iaga-code]_[date-time]_[publication-level].cdf
 
         Parameters:
         - observatory: IAGA code of the observatory.
         - date: Start date for the file.
         - type: Data type indicating the processing level.
         - interval: Data interval (e.g., 'minute', 'second').
-        - channels: List of channels (optional).
+        - channels: List of channels (optional). Not Implemented
 
         Returns:
         - The formatted file URL or path.
@@ -653,7 +670,7 @@ class ImagCDFFactory(TimeseriesFactory):
         # Get the publication level for the type
         publication_level = IMCDFPublicationLevel(data_type=type).to_string()
 
-        # Determine filename date format based on interval
+        # Format of Date/Time Portion of Filename based on interval see reference: https://tech-man.intermagnet.org/latest/appendices/dataformats.html#example-data-file:~:text=Format%20of%20Date,%EF%83%81
         if interval == "year":
             date_format = date.strftime("%Y")
         elif interval == "month":
@@ -667,10 +684,9 @@ class ImagCDFFactory(TimeseriesFactory):
         elif interval == "second":
             date_format = date.strftime("%Y%m%d_%H%M%S")
         else:
-            raise ValueError(f"Unsupported interval: {interval}")
+            raise ValueError(f"Unsupported interval: {interval}") #tenhertz currently not supported
 
-        # Default filename following ImagCDF convention
-        # Filename format: [iaga-code]_[date-time]_[publication-level].cdf
+        # Filename following ImagCDF convention, see reference: https://tech-man.intermagnet.org/latest/appendices/dataformats.html#imagcdf-file-names
         filename = f"{observatory.lower()}_{date_format}_{publication_level}.cdf"
 
         # If the urlTemplate explicitly specifies 'stdout', return 'stdout'
@@ -696,7 +712,7 @@ class ImagCDFFactory(TimeseriesFactory):
         }
 
         # Attempt to use the template provided in urlTemplate
-        if "{" in self.urlTemplate and "}" in self.urlTemplate:
+        if "{" in self.urlTemplate and "}" in self.urlTemplate: 
             try:
                 return self.urlTemplate.format(**params)
             except KeyError as e:
@@ -714,44 +730,6 @@ class ImagCDFFactory(TimeseriesFactory):
             f"Unsupported URL scheme in urlTemplate: {self.urlTemplate}"
         )
 
-    # Placeholder methods for interval and type names/abbreviations
-    def _get_interval_abbreviation(self, interval: DataInterval) -> str:
-        """Get the abbreviation for the data interval."""
-        abbreviations = {
-            "year": "yr",
-            "month": "mon",
-            "day": "day",
-            "hour": "hr",
-            "minute": "min",
-            "second": "sec",
-        }
-        return abbreviations.get(interval, "min")
-
-    def _get_interval_name(self, interval: DataInterval) -> str:
-        """Get the full name for the data interval."""
-        names = {
-            "year": "yearly",
-            "month": "monthly",
-            "day": "daily",
-            "hour": "hourly",
-            "minute": "minute",
-            "second": "second",
-        }
-        return names.get(interval, "minute")
-
-    def _get_type_name(self, type: DataType) -> str:
-        """Get the full name for the data type."""
-        type_names = {
-            "variation": "variation",
-            "definitive": "definitive",
-            "quasi-definitive": "quasi-definitive",
-            "provisional": "provisional",
-            "adjusted": "adjusted",
-            "none": "none",
-        }
-        return type_names.get(type, "variation")
-
-
     def _get_vector_elements(self):
         return {'X', 'Y', 'Z', 'H', 'D', 'E', 'V', 'I', 'F'}
     
diff --git a/geomagio/api/ws/Element.py b/geomagio/api/ws/Element.py
index ae5c4ead..25b876c5 100644
--- a/geomagio/api/ws/Element.py
+++ b/geomagio/api/ws/Element.py
@@ -45,3 +45,6 @@ ELEMENTS = [
 ]
 
 ELEMENT_INDEX = {e.id: e for e in ELEMENTS}
+
+
+TEMPERATURE_ELEMENTS_ID = [element.id for element in ELEMENTS if "Temperature" in element.name]
\ No newline at end of file
-- 
GitLab