Skip to content
Snippets Groups Projects
Commit 8b9a79fc authored by Shavers, Nicholas H's avatar Shavers, Nicholas H
Browse files

arcmin to degrees per cf convention. other minor updates per cf convention.

parent ff57e71f
No related branches found
No related tags found
1 merge request!379Netcdf mvp
...@@ -21,8 +21,8 @@ import numpy ...@@ -21,8 +21,8 @@ import numpy
M2R = numpy.pi / 180 / 60 # Minutes to Radians M2R = numpy.pi / 180 / 60 # Minutes to Radians
R2M = 180.0 / numpy.pi * 60 # Radians to Minutes 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.... # get geographic coordinates from....
# ### # ###
...@@ -434,3 +434,37 @@ def get_minutes_from_radians(r): ...@@ -434,3 +434,37 @@ def get_minutes_from_radians(r):
the radian value to be converted the radian value to be converted
""" """
return numpy.multiply(r, R2M) 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)
...@@ -6,6 +6,7 @@ import tempfile ...@@ -6,6 +6,7 @@ import tempfile
import shutil import shutil
import os import os
from geomagio import ChannelConverter
from geomagio.TimeseriesFactory import TimeseriesFactory from geomagio.TimeseriesFactory import TimeseriesFactory
from geomagio.api.ws.Element import ELEMENT_INDEX from geomagio.api.ws.Element import ELEMENT_INDEX
...@@ -142,8 +143,6 @@ class NetCDFFactory(TimeseriesFactory): ...@@ -142,8 +143,6 @@ class NetCDFFactory(TimeseriesFactory):
for ch in channels: for ch in channels:
channel_data[ch] = nc_dataset.variables[ch][:] channel_data[ch] = nc_dataset.variables[ch][:]
nc_dataset.close()
# If no time points, return empty # If no time points, return empty
if len(times_epoch) == 0: if len(times_epoch) == 0:
print("No time points found in the NetCDF data.") print("No time points found in the NetCDF data.")
...@@ -177,8 +176,19 @@ class NetCDFFactory(TimeseriesFactory): ...@@ -177,8 +176,19 @@ class NetCDFFactory(TimeseriesFactory):
# Create Trace # Create Trace
trace = Trace(data=np.array(data_array, dtype=float), header=stats) 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) stream.append(trace)
nc_dataset.close()
return stream return stream
def write_file_from_buffer(self, fh, timeseries: Stream, channels=None): def write_file_from_buffer(self, fh, timeseries: Stream, channels=None):
...@@ -198,6 +208,7 @@ class NetCDFFactory(TimeseriesFactory): ...@@ -198,6 +208,7 @@ class NetCDFFactory(TimeseriesFactory):
# Create an empty NetCDF structure with minimal CF compliance # Create an empty NetCDF structure with minimal CF compliance
with netCDF4.Dataset(fh, "w", format="NETCDF4") as nc_dataset: with netCDF4.Dataset(fh, "w", format="NETCDF4") as nc_dataset:
nc_dataset.title = "Geomagnetic Time Series Data" nc_dataset.title = "Geomagnetic Time Series Data"
nc_dataset.featureType = "timeSeries" # per cf conventions
nc_dataset.Conventions = "CF-1.6" nc_dataset.Conventions = "CF-1.6"
return return
...@@ -214,6 +225,7 @@ class NetCDFFactory(TimeseriesFactory): ...@@ -214,6 +225,7 @@ class NetCDFFactory(TimeseriesFactory):
print("No matching channels found after filtering.") print("No matching channels found after filtering.")
with netCDF4.Dataset(fh, "w", format="NETCDF4") as nc_dataset: with netCDF4.Dataset(fh, "w", format="NETCDF4") as nc_dataset:
nc_dataset.title = "Geomagnetic Time Series Data" nc_dataset.title = "Geomagnetic Time Series Data"
nc_dataset.featureType = "timeSeries" # per cf conventions
nc_dataset.Conventions = "CF-1.6" nc_dataset.Conventions = "CF-1.6"
return return
...@@ -239,6 +251,7 @@ class NetCDFFactory(TimeseriesFactory): ...@@ -239,6 +251,7 @@ class NetCDFFactory(TimeseriesFactory):
# Define global attributes with CF compliance # Define global attributes with CF compliance
nc_dataset.Conventions = "CF-1.6" nc_dataset.Conventions = "CF-1.6"
nc_dataset.title = "Geomagnetic Time Series Data" nc_dataset.title = "Geomagnetic Time Series Data"
nc_dataset.featureType = "timeSeries" # per cf conventions
nc_dataset.institution = getattr(stats, "agency_name", "") nc_dataset.institution = getattr(stats, "agency_name", "")
# nc_dataset.source = getattr(stats, "source", "") # nc_dataset.source = getattr(stats, "source", "")
history_entry = f"Created on {datetime.now(timezone.utc).isoformat()}" history_entry = f"Created on {datetime.now(timezone.utc).isoformat()}"
...@@ -268,6 +281,7 @@ class NetCDFFactory(TimeseriesFactory): ...@@ -268,6 +281,7 @@ class NetCDFFactory(TimeseriesFactory):
time_var.calendar = getattr(stats, "calendar", "standard") time_var.calendar = getattr(stats, "calendar", "standard")
time_var.standard_name = "time" time_var.standard_name = "time"
time_var.long_name = "Time" time_var.long_name = "Time"
time_var.axis = "T"
time_var[:] = times_epoch time_var[:] = times_epoch
# Assign valid range if available # Assign valid range if available
...@@ -278,6 +292,10 @@ class NetCDFFactory(TimeseriesFactory): ...@@ -278,6 +292,10 @@ class NetCDFFactory(TimeseriesFactory):
# Create channel variables with CF attributes # Create channel variables with CF attributes
for trace in timeseries: for trace in timeseries:
ch_name = trace.stats.channel.upper() 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: if len(ch_name) > 63:
raise ValueError( raise ValueError(
f"Channel name '{ch_name}' exceeds NetCDF variable name length limit." f"Channel name '{ch_name}' exceeds NetCDF variable name length limit."
...@@ -300,9 +318,13 @@ class NetCDFFactory(TimeseriesFactory): ...@@ -300,9 +318,13 @@ class NetCDFFactory(TimeseriesFactory):
# Set CF attributes # Set CF attributes
var.units = units var.units = units
var.long_name = getattr( if ELEMENT_INDEX.get(ch_name):
trace.stats, "long_name", ELEMENT_INDEX.get(ch_name).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( var.standard_name = getattr(
trace.stats, "standard_name", f"geomagnetic_field_{ch_name.lower()}" trace.stats, "standard_name", f"geomagnetic_field_{ch_name.lower()}"
) )
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment