feat: save start and end date of files

doc/source/notebooks/GRACE-Spatial-Error.ipynb

@@ -306,9 +306,9 @@
     "mm = np.arange(0, MMAX + 1)\n",
     "PLM2 = PLM[:, mm, :] ** 2\n",
     "# Calculating cos(m*phi)^2 and sin(m*phi)^2\n",
-    "phi = np.radians(grid.lon[np.newaxis, :])\n",
-    "ccos = np.cos(np.dot(mm[:, np.newaxis], phi)) ** 2\n",
-    "ssin = np.sin(np.dot(mm[:, np.newaxis], phi)) ** 2\n",
+    "phi = np.radians(grid.lon)\n",
+    "mp = np.einsum('m...,p...->mp...', mm, phi)\n",
+    "m_phi2 = (1.0 - 1j) * (np.exp(-2j * mp) + np.exp(2j * mp) + 2j) / 4.0\n",
     "\n",
     "# read load love numbers file\n",
     "# PREM outputs from Han and Wahr (1995)\n",
@@ -346,7 +346,7 @@
     "\n",
     "# dfactor is the degree dependent coefficients\n",
     "# for converting to spherical harmonic output units\n",
-    "factors = gravtk.units(lmax=LMAX).harmonic(hl, kl, ll).mmwe\n",
+    "factors = gravtk.units(lmax=LMAX).harmonic(hl, kl, ll)\n",
     "# mmwe, millimeters water equivalent\n",
     "dfactor = factors.get('mmwe')\n",
     "# units strings for output plots\n",
@@ -388,17 +388,14 @@
     "delta_Ylms = delta_Ylms.convolve(dfactor * wt)\n",
     "# smooth harmonics and convert to output units\n",
     "YLM2 = delta_Ylms.power(2.0).scale(1.0 / nsmth)\n",
+    "\n",
     "# Calculate fourier coefficients\n",
-    "d_cos = np.zeros((MMAX + 1, nlat))  # [m,th]\n",
-    "d_sin = np.zeros((MMAX + 1, nlat))  # [m,th]\n",
-    "# Calculating delta spatial values\n",
-    "for k in range(0, nlat):\n",
-    "    # summation over all spherical harmonic degrees\n",
-    "    d_cos[:, k] = np.sum(PLM2[:, :, k] * YLM2.clm, axis=0)\n",
-    "    d_sin[:, k] = np.sum(PLM2[:, :, k] * YLM2.slm, axis=0)\n",
+    "# summation over all spherical harmonic degrees\n",
+    "pconv2 = np.einsum('lmh...,lm...->mh...', PLM2, YLM2.ilm)\n",
     "\n",
-    "# Multiplying by c/s(phi#m) to get spatial maps (lon,lat)\n",
-    "grid.data = np.sqrt(np.dot(ccos.T, d_cos) + np.dot(ssin.T, d_sin)).T\n",
+    "# Multiplying by c/s(phi#m) to get spatial error map\n",
+    "# take the square root and drop imaginary component\n",
+    "grid.data = np.sqrt(np.einsum('mp...,mh...->hp...', m_phi2, pconv2)).real\n",
     "grid.mask = np.zeros_like(grid.data, dtype=bool)"
    ]
   },
@@ -514,11 +511,9 @@
   }
  ],
  "metadata": {
-  "interpreter": {
-   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
-  },
   "kernelspec": {
-   "display_name": "Python 3.8.10 64-bit",
+   "display_name": "py13",
+   "language": "python",
    "name": "python3"
   },
   "language_info": {
@@ -531,7 +526,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.6"
+   "version": "3.13.0"
   }
  },
  "nbformat": 4,

gravity_toolkit/grace_date.py

@@ -47,6 +47,7 @@
 
 UPDATE HISTORY:
     Updated 08/2026: added verbosity argument to print debugging information
+        use numpy datetime arrays to save start and end dates
     Updated 05/2023: use pathlib to define and operate on paths
     Updated 03/2023: use f-strings for formatting output date lines
         added regex formatting for CNES GRGS harmonics
@@ -110,7 +111,9 @@
 import pathlib
 import argparse
 import numpy as np
+import collections
 import gravity_toolkit as gravtk
+from datetime import datetime
 
 
 # PURPOSE: parses GRACE/GRACE-FO data files and assigns month numbers
@@ -177,85 +180,56 @@ def grace_date(base_dir, PROC='', DREL='', DSET='', OUTPUT=True, MODE=0o775):
     n_files = len(input_files)
 
     # define date variables
-    start_yr = np.zeros((n_files))  # year start date
-    end_yr = np.zeros((n_files))  # year end date
-    start_day = np.zeros((n_files))  # day number start date
-    end_day = np.zeros((n_files))  # day number end date
+    start_date = np.zeros((n_files), dtype='datetime64[s]')  # start date
+    end_date = np.zeros((n_files), dtype='datetime64[s]')  # end date
     mid_day = np.zeros((n_files))  # mid-month day
     tot_days = np.zeros((n_files))  # number of days since Jan 2002
     tdec = np.zeros((n_files))  # date in decimal form
     mon = np.zeros((n_files,), dtype=np.int64)  # GRACE/GRACE-FO month number
 
     # for each data file
     for t, infile in enumerate(input_files):
-        if PROC in (
-            'GRAZ',
-            'Swarm',
-        ):
+        if PROC in ('GRAZ', 'Swarm'):
             # get date lists for the start and end of fields
-            start_date, end_date = gravtk.time.parse_gfc_file(
+            start_date[t], end_date[t] = gravtk.time.parse_gfc_file(
                 infile, PROC, DSET
             )
-            # start and end year
-            start_yr[t] = np.float64(start_date[0])
-            end_yr[t] = np.float64(end_date[0])
-            # number of days in each month for the calendar year
-            dpm = gravtk.time.calendar_days(start_yr[t])
-            # start and end day of the year
-            start_day[t] = (
-                np.sum(dpm[: start_date[1] - 1])
-                + start_date[2]
-                + start_date[3] / 24.0
-                + start_date[4] / 1440.0
-                + start_date[5] / 86400.0
-            )
-            end_day[t] = (
-                np.sum(dpm[: end_date[1] - 1])
-                + end_date[2]
-                + end_date[3] / 24.0
-                + end_date[4] / 1440.0
-                + end_date[5] / 86400.0
-            )
         else:
             # get date lists for the start and end of fields
-            start_date, end_date = gravtk.time.parse_grace_file(infile)
-            # start and end year
-            start_yr[t] = np.float64(start_date[0])
-            end_yr[t] = np.float64(end_date[0])
-            # start and end day of the year
-            start_day[t] = np.float64(start_date[1])
-            end_day[t] = np.float64(end_date[1])
+            start_date[t], end_date[t] = gravtk.time.parse_grace_file(infile)
+        # start year and day of the year
+        start_struct = start_date[t].astype(datetime).timetuple()
+        start_yr = start_struct.tm_year
+        start_day = start_struct.tm_yday
+        # end year and day of the year
+        end_struct = end_date[t].astype(datetime).timetuple()
+        end_yr = end_struct.tm_year
+        end_day = end_struct.tm_yday
 
         # number of days in the starting year for leap and standard years
-        dpy = gravtk.time.calendar_days(start_yr[t]).sum()
+        dpy = gravtk.time.calendar_days(start_yr).sum()
         # end date taking into account measurements taken on different years
-        end_cyclic = (end_yr[t] - start_yr[t]) * dpy + end_day[t]
+        end_cyclic = (end_yr - start_yr) * dpy + end_day
         # calculate mid-month value
-        mid_day[t] = np.mean([start_day[t], end_cyclic])
+        mid_day[t] = np.mean([start_day, end_cyclic])
 
         # calculate Modified Julian Day from start_yr and mid_day
         MJD = gravtk.time.convert_calendar_dates(
-            start_yr[t], 1.0, mid_day[t], epoch=(1858, 11, 17, 0, 0, 0)
+            start_yr, 1.0, mid_day[t], epoch=(1858, 11, 17, 0, 0, 0)
         )
         # convert from Modified Julian Days to calendar dates
         cal_date = gravtk.time.convert_julian(MJD + 2400000.5)
 
         # Calculating the mid-month date in decimal form
-        tdec[t] = start_yr[t] + mid_day[t] / dpy
+        tdec[t] = start_yr + mid_day[t] / dpy
 
         # Calculation of total days since start of campaign
-        count = 0
-        n_yrs = np.int64(start_yr[t] - 2002)
-        # for each of the GRACE years up to the file year
-        for iyr in range(n_yrs):
-            # year
-            year = 2002 + iyr
-            # add all days from prior years to count
-            # number of days in year i (if leap year or standard year)
-            count += gravtk.time.calendar_days(year).sum()
-
-        # calculating the total number of days since 2002
-        tot_days[t] = np.mean([count + start_day[t], count + end_cyclic])
+        # add all days from prior years to count
+        # number of days in year i (if leap year or standard year)
+        count = np.sum(
+            [gravtk.time.calendar_days(y) for y in range(2002, start_yr)]
+        )
+        tot_days[t] = count + mid_day[t]
 
         # Calculates the month number (or 10-day number for CNES RL01,RL02)
         if (PROC == 'CNES') and (DREL in ('RL01', 'RL02')):
@@ -285,18 +259,33 @@ def grace_date(base_dir, PROC='', DREL='', DSET='', OUTPUT=True, MODE=0o775):
         print('{0} {1:>10} {2:>11} {3:>10} {4:>13}'.format(*args), file=fid)
 
     # create python dictionary mapping input file names with GRACE months
-    grace_files = {}
+    grace_files = collections.OrderedDict()
+    # add attributes
+    grace_files.attrs = {
+        'start_date': start_date,
+        'end_date': end_date,
+        'time_decimal': tdec,
+    }
     # for each data file
     for t, infile in enumerate(input_files):
         # add file to python dictionary mapped to GRACE/GRACE-FO month
         grace_files[mon[t]] = grace_dir.joinpath(infile)
         # print to GRACE dates ascii file (NOTE: tot_days will be rounded)
         if OUTPUT:
+            # start year and day of the year
+            start_struct = start_date[t].astype(datetime).timetuple()
+            start_yr = start_struct.tm_year
+            start_day = start_struct.tm_yday
+            # end year and day of the year
+            end_struct = end_date[t].astype(datetime).timetuple()
+            end_yr = end_struct.tm_year
+            end_day = end_struct.tm_yday
+            # print to GRACE dates ascii file
             print(
                 (
                     f'{tdec[t]:13.8f} {mon[t]:03d} '
-                    f'{start_yr[t]:8.0f} {start_day[t]:03.0f} '
-                    f'{end_yr[t]:8.0f} {end_day[t]:03.0f} '
+                    f'{start_yr:8.0f} {start_day:03.0f} '
+                    f'{end_yr:8.0f} {end_day:03.0f} '
                     f'{tot_days[t]:8.0f}'
                 ),
                 file=fid,

gravity_toolkit/grace_input_months.py

@@ -110,6 +110,7 @@
 
 UPDATE HISTORY:
     Updated 08/2026: use upstream file logger for verbose output
+        save start and end date of files as numpy datetime64[s] arrays
     Updated 10/2023: standardize ocean model for UCI degree 1 coefficients
     Updated 05/2023: use pathlib to define and operate on paths
     Updated 04/2023: use release-03 GFZ GravIS SLR and geocenter files
@@ -363,7 +364,7 @@ def grace_input_months(
     n_cons = len(months)
 
     # Initializing input data matrices
-    grace_Ylms = {}
+    grace_Ylms = collections.OrderedDict()
     grace_Ylms['clm'] = np.zeros((LMAX + 1, MMAX + 1, n_cons))
     grace_Ylms['slm'] = np.zeros((LMAX + 1, MMAX + 1, n_cons))
     grace_Ylms['eclm'] = np.zeros((LMAX + 1, MMAX + 1, n_cons))
@@ -373,6 +374,11 @@ def grace_input_months(
     # output dimensions
     grace_Ylms['l'] = np.arange(LMAX + 1)
     grace_Ylms['m'] = np.arange(MMAX + 1)
+    # datetime arrays for start and end dates
+    grace_Ylms.attrs = {
+        'start_date': np.zeros((n_cons), dtype='datetime64[s]'),
+        'end_date': np.zeros((n_cons), dtype='datetime64[s]'),
+    }
 
     # attributes for processing run
     attributes = collections.OrderedDict()
@@ -383,11 +389,14 @@ def grace_input_months(
     grace_files = grace_date(
         base_dir, PROC=PROC, DREL=DREL, DSET=DSET, OUTPUT=False
     )
+    # get the list of GRACE/GRACE-FO months for the input files
+    file_months = list(grace_files.keys())
 
     # importing data from GRACE/GRACE-FO files
     for i, grace_month in enumerate(months):
         # read spherical harmonic data products
         infile = grace_files[grace_month]
+        j = file_months.index(grace_month)
         # log input file if debugging
         logger.debug(f'Reading file {i:d}: {str(infile)}')
         # read GRACE/GRACE-FO/Swarm file
@@ -408,6 +417,9 @@ def grace_input_months(
         # copy date variables
         grace_Ylms['time'][i] = np.copy(Ylms['time'])
         grace_Ylms['month'][i] = np.int64(grace_month)
+        # copy start and end dates from grace_files attributes
+        grace_Ylms.attrs['start_date'][i] = grace_files.attrs['start_date'][j]
+        grace_Ylms.attrs['end_date'][i] = grace_files.attrs['end_date'][j]
         # copy input file basename
         attributes['lineage'][i] = pathlib.Path(infile).stem