From 8b9a79fc3ea7cf84595e95d68512779c5da2c425 Mon Sep 17 00:00:00 2001
From: Nicholas Shavers <nshavers@contractor.usgs.gov>
Date: Mon, 6 Jan 2025 10:52:26 -0800
Subject: [PATCH] arcmin to degrees per cf convention. other minor updates per
cf convention.
---
geomagio/ChannelConverter.py | 38 ++++++++++++++++++++++++++++++--
geomagio/netcdf/NetCDFFactory.py | 32 ++++++++++++++++++++++-----
2 files changed, 63 insertions(+), 7 deletions(-)
diff --git a/geomagio/ChannelConverter.py b/geomagio/ChannelConverter.py
index 0ca91bfd..3fc705ff 100644
--- a/geomagio/ChannelConverter.py
+++ b/geomagio/ChannelConverter.py
@@ -21,8 +21,8 @@ import numpy
M2R = numpy.pi / 180 / 60 # Minutes to Radians
R2M = 180.0 / numpy.pi * 60 # Radians to Minutes
-
-
+M2D = 1.0 / 60.0 # Minutes to Degrees
+D2M = 60.0 # Degrees to Minutes
# ###
# get geographic coordinates from....
# ###
@@ -434,3 +434,37 @@ def get_minutes_from_radians(r):
the radian value to be converted
"""
return numpy.multiply(r, R2M)
+
+
+def get_degrees_from_minutes(m):
+ """
+ Converts angular minutes to degrees.
+
+ Parameters
+ ----------
+ m : float or array_like
+ The angular minutes to be converted.
+
+ Returns
+ -------
+ float or array_like
+ The corresponding degrees.
+ """
+ return numpy.multiply(m, M2D)
+
+
+def get_minutes_from_degrees(d):
+ """
+ Converts degrees to angular minutes.
+
+ Parameters
+ ----------
+ d : float or array_like
+ The degrees to be converted.
+
+ Returns
+ -------
+ float or array_like
+ The corresponding angular minutes.
+ """
+ return numpy.multiply(d, D2M)
diff --git a/geomagio/netcdf/NetCDFFactory.py b/geomagio/netcdf/NetCDFFactory.py
index 590e1ead..3a92d580 100644
--- a/geomagio/netcdf/NetCDFFactory.py
+++ b/geomagio/netcdf/NetCDFFactory.py
@@ -6,6 +6,7 @@ import tempfile
import shutil
import os
+from geomagio import ChannelConverter
from geomagio.TimeseriesFactory import TimeseriesFactory
from geomagio.api.ws.Element import ELEMENT_INDEX
@@ -142,8 +143,6 @@ class NetCDFFactory(TimeseriesFactory):
for ch in channels:
channel_data[ch] = nc_dataset.variables[ch][:]
- nc_dataset.close()
-
# If no time points, return empty
if len(times_epoch) == 0:
print("No time points found in the NetCDF data.")
@@ -177,8 +176,19 @@ class NetCDFFactory(TimeseriesFactory):
# Create Trace
trace = Trace(data=np.array(data_array, dtype=float), header=stats)
+
+ # Convert minutes to arcradians
+ if channel_name in ["D", "I"]:
+ if nc_dataset.variables[channel_name].units.startswith("deg"):
+ trace.data = ChannelConverter.get_minutes_from_degrees(trace.data)
+ elif nc_dataset.variables[channel_name].units.startswith("rad"):
+ trace.data = ChannelConverter.get_minutes_from_radians(trace.data)
+ elif nc_dataset.variables[channel_name].units.startswith("arc"):
+ pass # nothing to do.
+
stream.append(trace)
+ nc_dataset.close()
return stream
def write_file_from_buffer(self, fh, timeseries: Stream, channels=None):
@@ -198,6 +208,7 @@ class NetCDFFactory(TimeseriesFactory):
# Create an empty NetCDF structure with minimal CF compliance
with netCDF4.Dataset(fh, "w", format="NETCDF4") as nc_dataset:
nc_dataset.title = "Geomagnetic Time Series Data"
+ nc_dataset.featureType = "timeSeries" # per cf conventions
nc_dataset.Conventions = "CF-1.6"
return
@@ -214,6 +225,7 @@ class NetCDFFactory(TimeseriesFactory):
print("No matching channels found after filtering.")
with netCDF4.Dataset(fh, "w", format="NETCDF4") as nc_dataset:
nc_dataset.title = "Geomagnetic Time Series Data"
+ nc_dataset.featureType = "timeSeries" # per cf conventions
nc_dataset.Conventions = "CF-1.6"
return
@@ -239,6 +251,7 @@ class NetCDFFactory(TimeseriesFactory):
# Define global attributes with CF compliance
nc_dataset.Conventions = "CF-1.6"
nc_dataset.title = "Geomagnetic Time Series Data"
+ nc_dataset.featureType = "timeSeries" # per cf conventions
nc_dataset.institution = getattr(stats, "agency_name", "")
# nc_dataset.source = getattr(stats, "source", "")
history_entry = f"Created on {datetime.now(timezone.utc).isoformat()}"
@@ -268,6 +281,7 @@ class NetCDFFactory(TimeseriesFactory):
time_var.calendar = getattr(stats, "calendar", "standard")
time_var.standard_name = "time"
time_var.long_name = "Time"
+ time_var.axis = "T"
time_var[:] = times_epoch
# Assign valid range if available
@@ -278,6 +292,10 @@ class NetCDFFactory(TimeseriesFactory):
# Create channel variables with CF attributes
for trace in timeseries:
ch_name = trace.stats.channel.upper()
+ # Convert arcminutes to radians
+ if "D" == ch_name:
+ trace.data = ChannelConverter.get_degrees_from_minutes(trace.data)
+
if len(ch_name) > 63:
raise ValueError(
f"Channel name '{ch_name}' exceeds NetCDF variable name length limit."
@@ -300,9 +318,13 @@ class NetCDFFactory(TimeseriesFactory):
# Set CF attributes
var.units = units
- var.long_name = getattr(
- trace.stats, "long_name", ELEMENT_INDEX.get(ch_name).name
- )
+ if ELEMENT_INDEX.get(ch_name):
+ # remove arcminutes, other things
+ long_name = ELEMENT_INDEX.get(ch_name).name.replace(
+ "(arcminutes)", ""
+ )
+ long_name = long_name.strip()
+ var.long_name = getattr(trace.stats, "long_name", long_name)
var.standard_name = getattr(
trace.stats, "standard_name", f"geomagnetic_field_{ch_name.lower()}"
)
--
GitLab