From a7fdf34f0ed5d445efc1baab3a182a717a503f61 Mon Sep 17 00:00:00 2001 From: Peter Norris Date: Wed, 3 Apr 2024 12:46:09 -0400 Subject: [PATCH 1/4] pmn: initial commit --- .../NRLSSI2/Mg_SB_from_daily_file.py | 159 +++++++++++ GEOS_RadiationShared/NRLSSI2/README.md | 14 + .../NRLSSI2/TSI_Mg_SB_merged_from_daily.py | 158 +++++++++++ .../NRLSSI2/TSI_from_daily_files.py | 249 ++++++++++++++++++ 4 files changed, 580 insertions(+) create mode 100644 GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py create mode 100644 GEOS_RadiationShared/NRLSSI2/README.md create mode 100644 GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py create mode 100644 GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py diff --git a/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py b/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py new file mode 100644 index 0000000..727b77f --- /dev/null +++ b/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py @@ -0,0 +1,159 @@ +''' +Get Solar Mg and SB indices from NRLSSI2 daily auxilliary text file +IMPORTANT: The filenames list below will be read in order, and any date that +already has final data will not be overwriten by data for that date in a later +file. For this reason, the filenames list should ONLY be APPENDED to, which +will ensure that any older final data is not overwritten, thus ensuring +historical reproducability. +''' + +import os +import numpy as np +import matplotlib.pyplot as plt +from calendar import monthrange +from datetime import datetime, timedelta + +# because datetime.strftime doesnt work before 1900 +def _yyyymmdd(t): + return('%04d%02d%02d' % (t.year, t.month, t.day)) + +class Mg_SB_Daily: + ''' Daily time series for Solar Mg and SB indices from NRLSSI2 ''' + + def __init__(self, + filenames=[ + 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20170630_c20170928.txt', + 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20180331_c20180418.txt', + 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20191231_c20200225.txt', + 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20201231_c20210204.txt', + 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20220630_c20220803.txt', + 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20230331_c20230411.txt', + 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20231231_c20240221.txt', + ], + verbose=False): + + self.filenames = filenames + self.verbose = verbose + + NB = os.environ['NOBACKUP'] + + if verbose: + print('reading Mg and SB indicies from', filenames) + + # read the files + self._data = {} + self.nfinal = 0 + first = True + for filename in filenames: + with open(os.sep.join((NB,'NRLSSI2',filename))) as f: + lines = f.readlines() + + # process the data, ignoring the header (top line) + for line in lines[1:]: + datestr, SB, Mg, status = line.strip().split(',') + # yyyy-mm-dd to yyyymmdd + yyyymmdd = datestr.replace('-','') + # make data real + SB = float(SB) + Mg = float(Mg) + # validate status to rc + if status == 'final': + rc = 0 + elif status == 'prelim': + rc = -1 + else: + raise RuntimeError('invalid status detected: %s' % status) + # if a date is already finalized, ignore further records + # for that date to ensure historical reproducability + if yyyymmdd in self._data: + rc_existing = self._data[yyyymmdd][0] + if rc_existing == 0: continue + # load data dictionary + self._data[yyyymmdd] = (rc, Mg, SB) + # keep track of min and max final times + if rc == 0: + self.nfinal += 1 + d = datetime.strptime(yyyymmdd,'%Y%m%d').date() + if first: + self.date_min_final = d + self.date_max_final = d + first = False + else: + if d < self.date_min_final: self.date_min_final = d + if d > self.date_max_final: self.date_max_final = d + + def viewkeys(self): + return self._data.viewkeys() + + def getday(self, yyyymmdd): + ''' + get daily Mg and SB values for daily string yyyymmdd + returns (rc, Mg, SB) + rc = 0 (final), -1 (prelim), -2 (unavailable) + ''' + if yyyymmdd not in self._data: + return -2, np.nan, np.nan + else: + return self._data[yyyymmdd] + + def gettime(self, t): + ''' + get time-interpolated Mg and SB values for datetime t + returns (rc, Mg, SB) + rc = 0 (final), -1 (prelim), -2 (unavailable) + ''' + + # assume daily average valid at noon + tnoon = datetime(t.year,t.month,t.day,12) + vnoon = self.getday(_yyyymmdd(tnoon)) + if t == tnoon: return vnoon + + # other noon bracketing t + tother = tnoon + timedelta(days=(-1 if t < tnoon else +1)) + vother = self.getday(_yyyymmdd(tother)) + + # fraction that the other daily average contributes + fother = abs((t-tnoon).total_seconds()) / 86400. + + # only interpolate if both days exist + if vnoon[0] == -2 or vother[0] == -2: + return (-2, np.nan, np.nan) + else: + # now both days available + Mg = vnoon[1] * (1-fother) + vother[1] * fother + SB = vnoon[2] * (1-fother) + vother[2] * fother + rc = -1 if vnoon[0] == -1 or vother[0] == -1 else 0 + return (rc, Mg, SB) + + def _plot_all_final(self): + spyear = 365.25 * 86400. + d0 = self.date_min_final; d = d0 + dy = []; Mg = []; SB = [] + while (d <= self.date_max_final): + v = self.getday(_yyyymmdd(d)) + if v[0] == 0: + # only plot finals + dy.append((d-d0).total_seconds()/spyear) + Mg.append(v[1]) + SB.append(v[2]) + d += timedelta(days=1) + plt.subplot(211); plt.plot(dy,Mg,'b-'); plt.ylabel('Mg') + plt.subplot(212); plt.plot(dy,SB,'b-'); plt.ylabel('SB') + plt.xlabel('years since %s' % self.date_min_final) + +if __name__ == '__main__': + + MgSB = Mg_SB_Daily() + print('16000101', MgSB.getday('16000101')) + print('20161130', MgSB.getday('20161130')) + print('20161201', MgSB.getday('20161201')) + print('20161202', MgSB.getday('20161202')) + t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S') + print(t, MgSB.gettime(t)) + + plt.figure() + MgSB._plot_all_final() + plt.savefig(os.sep.join(('gx','MgSB_plot_all_final.png')), + pad_inches=0.33,bbox_inches='tight',dpi=100) + plt.show() + diff --git a/GEOS_RadiationShared/NRLSSI2/README.md b/GEOS_RadiationShared/NRLSSI2/README.md new file mode 100644 index 0000000..bc012c2 --- /dev/null +++ b/GEOS_RadiationShared/NRLSSI2/README.md @@ -0,0 +1,14 @@ +get data files from https://www.ncei.noaa.gov/data/total-solar-irradiance/access/ +Note: can use wget https://... from discover +1. from ancilliary-data/, e.g., tsi-ssi_v02r01_model-input-time-series_s18820101_e20180331_c20180418.txt + !!! NOTE: Please read comments in docstring of Mg_SB_from_daily_file.py BEFORE downloading !!! + In particular, + -- the filenames list in that file should ONLY be APPENDED to, + which will ensure that any existing final data is not overwritten, + thus ensuring historical reproduceability. +2. from daily/, e.g., tsi_v02r01_daily_s20170101_e20171231_c20180227.nc + !!! NOTE: Please read comments in docstring of TSI_from_daily_files.py BEFORE downloading !!! + In particular, + -- there should be no files with overlapping time periods, and + -- you should NEVER overwrite any existing non-preliminary file, + since this may cause historical non-reproduceability. diff --git a/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py b/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py new file mode 100644 index 0000000..b0c6667 --- /dev/null +++ b/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py @@ -0,0 +1,158 @@ +''' +Get TSI and Solar Mg and SB indices from NRLSSI2 daily data +''' + +import os +import numpy as np +import matplotlib.pyplot as plt +from calendar import monthrange +from datetime import datetime, timedelta +from TSI_from_daily_files import TSI_Daily +from Mg_SB_from_daily_file import Mg_SB_Daily + +# because datetime.strftime doesnt work before 1900 +def _yyyymmdd(t): + return('%04d%02d%02d' % (t.year, t.month, t.day)) + +class TSI_Mg_SB_Merged_Daily: + ''' Daily time series for TSI and Solar Mg and SB indices from NRLSSI2 ''' + + def __init__(self, verbose=False): + + self.verbose = verbose + + # get each data set + zM = Mg_SB_Daily(verbose=verbose) + zT = TSI_Daily(verbose=verbose) + + # form the INTERSECTION: + # (each set is indexed by yyyymmdd, unique in each set) + self._data = {} + self.nfinal = 0 + first = True + for yyyymmdd in (zM.viewkeys() & zT.viewkeys()): + + # load data dictionary + rc1, Mg, SB = zM.getday(yyyymmdd) + rc2, TSI, TSI_UNC = zT.getday(yyyymmdd) + rc = -1 if rc1 == -1 or rc2 == -1 else 0 + self._data[yyyymmdd] = (rc, TSI, Mg, SB, TSI_UNC) + + # keep track of min and max final times + if rc == 0: + self.nfinal += 1 + d = datetime.strptime(yyyymmdd,'%Y%m%d').date() + if first: + self.date_min_final = d + self.date_max_final = d + first = False + else: + if d < self.date_min_final: self.date_min_final = d + if d > self.date_max_final: self.date_max_final = d + + def getday(self, yyyymmdd): + ''' + get daily TSI, Mg and SB values for daily string yyyymmdd + returns (rc, TSI, Mg, SB, TSI_UNC) + rc = 0 (final), -1 (prelim), -2 (unavailable) + ''' + if yyyymmdd not in self._data: + return -2, np.nan, np.nan, np.nan + else: + return self._data[yyyymmdd] + + def gettime(self, t): + ''' + get time-interpolated TSI, Mg, and SB values for datetime t + returns (rc, TSI, Mg, SB, TSI_UNC) + rc = 0 (final), -1 (prelim), -2 (unavailable) + ''' + + # assume daily average valid at noon + tnoon = datetime(t.year,t.month,t.day,12) + vnoon = self.getday(_yyyymmdd(tnoon)) + if t == tnoon: return vnoon + + # other noon bracketing t + tother = tnoon + timedelta(days=(-1 if t < tnoon else +1)) + vother = self.getday(_yyyymmdd(tother)) + + # fraction that the other daily average contributes + fother = abs((t-tnoon).total_seconds()) / 86400. + + # only interpolate if both days exist + if vnoon[0] == -2 or vother[0] == -2: + return (-2, np.nan, np.nan, np.nan) + else: + # now both days available + TSI = vnoon[1] * (1-fother) + vother[1] * fother + Mg = vnoon[2] * (1-fother) + vother[2] * fother + SB = vnoon[3] * (1-fother) + vother[3] * fother + TSI_UNC = vnoon[4] * (1-fother) + vother[4] * fother + rc = -1 if vnoon[0] == -1 or vother[0] == -1 else 0 + return (rc, TSI, Mg, SB, TSI_UNC) + + def final_date_range(self, + yyyymmdd0, yyyymmdd1, + no_gaps=True): + ''' final only lists over date range ''' + + d0 = datetime.strptime(yyyymmdd0,'%Y%m%d').date() + d1 = datetime.strptime(yyyymmdd1,'%Y%m%d').date() + d = d0; dates = []; TSI = []; Mg = []; SB = [] + while (d <= d1): + v = self.getday(_yyyymmdd(d)) + if v[0] == 0: + # only include finals + dates.append(d) + TSI.append(v[1]) + Mg.append( v[2]) + SB.append( v[3]) + elif no_gaps: + raise RuntimeError('no_gaps: intervening non-final date found') + d += timedelta(days=1) + return dates, TSI, Mg, SB + + def _plot_all_final(self): + dates, TSI, Mg, SB = self.final_date_range( + _yyyymmdd(self.date_min_final), + _yyyymmdd(self.date_max_final), + no_gaps=False) + spyear = 365.25 * 86400. + dy = [(d-self.date_min_final).total_seconds()/spyear for d in dates] + plt.subplot(311); plt.plot(dy,TSI,'b-'); plt.ylabel('TSI') + plt.subplot(312); plt.plot(dy, Mg,'b-'); plt.ylabel('Mg') + plt.subplot(313); plt.plot(dy, SB,'b-'); plt.ylabel('SB') + plt.xlabel('years since %s' % self.date_min_final) + + def output_final_textfile(self, filename): + f = open(filename,'w') + f.write('# NRLSSI2 daily input\n') + f.write('# treat daily values as valid at 12:00 GMT\n') + f.write('# yyyy doy TSI:W/m2 MgIndex SBindex\n') + for d, TSI, Mg, SB in zip(*self.final_date_range( + _yyyymmdd(self.date_min_final), _yyyymmdd(self.date_max_final))): + f.write(' %04d %03d %8.3f %8.6f %9.4f\n' % ( + d.year, d.timetuple().tm_yday, TSI, Mg, SB)) + f.close() + +if __name__ == '__main__': + + NB = os.environ['NOBACKUP'] + + z = TSI_Mg_SB_Merged_Daily() + print('16000101', z.getday('16000101')) + print('20161130', z.getday('20161130')) + print('20161201', z.getday('20161201')) + print('20161202', z.getday('20161202')) + t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S') + print(t, z.gettime(t)) + + z.output_final_textfile(os.sep.join((NB,'NRLSSI2','output','NRLSSI2.txt'))) + + plt.figure(figsize=(6,12)) + z._plot_all_final() + plt.savefig(os.sep.join(('gx','TSInMgSB_plot_all_final.png')), + pad_inches=0.33,bbox_inches='tight',dpi=100) + plt.show() + diff --git a/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py b/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py new file mode 100644 index 0000000..2f47947 --- /dev/null +++ b/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py @@ -0,0 +1,249 @@ +''' +Get TSI from daily NRLSSI2 netCDF files +============================================================ +IMPORTANT RULES THAT GOVERN data/tsi_*_daily_*.nc FILES: +============================================================ +In order to maintain historical reproducabilty and not fail +the error checking below, please ensure the following for +data/tsi_*_daily_*.nc files: +1. files specified in the file_template argument must +obey 'data/tsi_*_daily_s????????_e????????_c????????.nc' +2. the first '*' is the version e.g., 'v02r01', with an +optional trailing '-preliminary' +3. The starting date s???????? should be of the form + sYYYYMMDD, where MMDD will typically be 0101. +4. The ending date e???????? should be of the form + eYYYYMMDD, where MMDD will typically be 1231, + but may be before then for preliminary files. +3. The creation date c???????? should be of the form + cYYYYMMDD. +5. It is up to the downloader and maintainer of these files +to ENSURE that + a. There should be no files with overlapping time periods + (For example, multiple versions / creation dates for + the same start and end date are not allowed. You must + keep only one, and probably the ORIGINAL(!) one, as + discussed in point (c.) below.) + b. You should NEVER overwrite an existing non-preliminary + file, since this may lose historical reproducability. + (Preliminary files may be overwritten with non-prelim- + inary ones, or with preliminary files of a later + creation date.) + c. Even if a new version comes along, that should only be + used for new years (or else you should maintain data + directories for different versions). Overwriting an + older version with a new one for an existing period in + this data directory, will destroy historical reproduc- + ability. +============================================================ +''' + +import os +import sys +import numpy as np +from glob import glob +from netCDF4 import Dataset +import matplotlib.pyplot as plt +from calendar import monthrange +from datetime import datetime, timedelta +from operator import itemgetter + +def _yyyymmdd(t): + ''' because datetime.strftime doesnt work before 1900 ''' + return('%04d%02d%02d' % (t.year, t.month, t.day)) + +class TSI_Daily: + ''' Daily time series for NRLSSI2 TSI ''' + + def __init__(self, + file_template='data/tsi_*_daily_s????????_e????????_c????????.nc', + verbose=False): + + self.file_template = file_template + self.verbose = verbose + + if verbose: + print('reading TSI from', file_template) + + # ============== + # Validate files + # ============== + + NB = os.environ['NOBACKUP'] + + # collect templated files and validate + files = [] + for path in glob(os.sep.join((NB,'NRLSSI2',file_template))): + filename = os.path.basename(path) + # validate filename + v = filename.split('_') + if v[0] != 'tsi' or v[2] != 'daily': + raise RuntimeError('%s not an NRLSSI2 tsi daily file' % file) + # version preliminary? + version = v[1] + if '-preliminary' in version: + rc = -1 + else: + rc = 0 + # extract inclusive date range of file + ds = datetime.strptime(v[3],'s%Y%m%d').date() + de = datetime.strptime(v[4],'e%Y%m%d').date() + if ds > de: + raise RuntimeError('bad timestamp ordering in %s' % filename) + # add to files to process + files.append((ds, de, rc, version, path)) + + if len(files) == 0: + raise RuntimeError('no TSI files to process') + + + # ======================== + # Fail if any time overlap + # ======================== + + # sort the files by start date + files.sort(key=itemgetter(0)) + + # detect any overlapping date ranges + if len(files) > 1: + ds0, de0, r0, v0, f0 = files[0] + for ds1, de1, r1, v1, f1 in files[1:]: + if ds0 == ds1: + # same start so overlapping + raise RuntimeError('date overlap: %s, %s' % (f0, f1)) + else: + # now ds0 < ds1 because of sort + if de0 >= ds1: + raise RuntimeError('date overlap: %s, %s' % (f0, f1)) + ds0, de0, r0, v0, f0 = ds1, de1, r1, v1, f1 + + # ================= + # Process the files + # ================= + + self._data = {} + self.nfinal = 0 + first = True + for ds, de, rc, version, path in files: + + # open file + if self.verbose: + print('processing path %s' % path) + hf = Dataset(path,'r') + + # extract dates + htime = hf.variables['time'] + tref = datetime.strptime(htime.units,'days since %Y-%m-%d %H:%M:%S') + dates = [(tref+timedelta(days=float(days))).date() for days in htime[:]] + + # read TSI, TSI_UNC + TSIs = hf.variables['TSI' ][:].astype(float) + TSI_UNCs = hf.variables['TSI_UNC'][:].astype(float) + + # close file + hf.close() + if not self.verbose: + sys.stdout.write('.') + sys.stdout.flush() + + # process file + for d, TSI, TSI_UNC in zip(dates, TSIs, TSI_UNCs): + + # load data dictionary + yyyymmdd = _yyyymmdd(d) + if yyyymmdd in self._data: + raise RuntimeError('date duplicates detected for %s' % yyyymmdd) + else: + self._data[yyyymmdd] = (rc, TSI, TSI_UNC) + + # keep track of min and max final times + if rc == 0: + self.nfinal += 1 + if first: + self.date_min_final = d + self.date_max_final = d + first = False + else: + if d < self.date_min_final: self.date_min_final = d + if d > self.date_max_final: self.date_max_final = d + + if not self.verbose: + sys.stdout.write('\n') + sys.stdout.flush() + + def viewkeys(self): + return self._data.viewkeys() + + def getday(self, yyyymmdd): + ''' + get daily TSI values for daily string yyyymmdd + returns (rc, TSI, TSI_UNC) + rc = 0 (final), -1 (prelim), -2 (unavailable) + ''' + if yyyymmdd not in self._data: + return -2, np.nan, np.nan + else: + return self._data[yyyymmdd] + + def gettime(self, t): + ''' + get time-interpolated TSI values for datetime t + returns (rc, TSI, TSI_UNC) + rc = 0 (final), -1 (prelim), -2 (unavailable) + ''' + + # assume daily average valid at noon + tnoon = datetime(t.year,t.month,t.day,12) + vnoon = self.getday(_yyyymmdd(tnoon)) + if t == tnoon: return vnoon + + # other noon bracketing t + tother = tnoon + timedelta(days=(-1 if t < tnoon else +1)) + vother = self.getday(_yyyymmdd(tother)) + + # fraction that the other daily average contributes + fother = abs((t-tnoon).total_seconds()) / 86400. + + # only interpolate if both days exist + if vnoon[0] == -2 or vother[0] == -2: + return (-2, np.nan, np.nan) + else: + # now both days available + TSI = vnoon[1] * (1-fother) + vother[1] * fother + TSI_UNC = vnoon[2] * (1-fother) + vother[2] * fother + rc = -1 if vnoon[0] == -1 or vother[0] == -1 else 0 + return (rc, TSI, TSI_UNC) + + def _plot_all_final(self): + spyear = 365.25 * 86400. + d0 = self.date_min_final; d = d0 + dy = []; TSI = []; TSI_UNC = [] + while (d <= self.date_max_final): + v = self.getday(_yyyymmdd(d)) + if v[0] == 0: + # only plot finals + dy.append((d-d0).total_seconds()/spyear) + TSI.append(v[1]) + TSI_UNC.append(v[2]) + d += timedelta(days=1) + plt.subplot(211); plt.plot(dy,TSI,'b-'); plt.ylabel('TSI [W/m2]') + plt.title(r'%s'%self.file_template,fontsize=10) + plt.subplot(212); plt.plot(dy,TSI_UNC,'b-'); plt.ylabel('TSI_UNC [W/m2]') + plt.xlabel('years since %s' % self.date_min_final) + +if __name__ == '__main__': + + TSI = TSI_Daily() + print('16000101', TSI.getday('16000101')) + print('20161130', TSI.getday('20161130')) + print('20161201', TSI.getday('20161201')) + print('20161202', TSI.getday('20161202')) + t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S') + print(t, TSI.gettime(t)) + + plt.figure() + TSI._plot_all_final() + plt.savefig(os.sep.join(('gx','TSI_plot_all_final.png')), + pad_inches=0.33,bbox_inches='tight',dpi=100) + plt.show() + From 0e6f3f77fe6b5a7dce2334ae1fe4b511ba8b12ae Mon Sep 17 00:00:00 2001 From: Peter Norris Date: Wed, 3 Apr 2024 14:14:29 -0400 Subject: [PATCH 2/4] pmn: python3 adaptation completed --- .../NRLSSI2/Mg_SB_from_daily_file.py | 20 +++++++++---------- .../NRLSSI2/TSI_Mg_SB_merged_from_daily.py | 18 ++++++++--------- .../NRLSSI2/TSI_from_daily_files.py | 20 +++++++++---------- 3 files changed, 29 insertions(+), 29 deletions(-) diff --git a/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py b/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py index 727b77f..e3735a3 100644 --- a/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py +++ b/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py @@ -82,8 +82,8 @@ def __init__(self, if d < self.date_min_final: self.date_min_final = d if d > self.date_max_final: self.date_max_final = d - def viewkeys(self): - return self._data.viewkeys() + def keys(self): + return self._data.keys() def getday(self, yyyymmdd): ''' @@ -144,16 +144,16 @@ def _plot_all_final(self): if __name__ == '__main__': MgSB = Mg_SB_Daily() - print('16000101', MgSB.getday('16000101')) - print('20161130', MgSB.getday('20161130')) - print('20161201', MgSB.getday('20161201')) - print('20161202', MgSB.getday('20161202')) - t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S') - print(t, MgSB.gettime(t)) +# print('16000101', MgSB.getday('16000101')) +# print('20161130', MgSB.getday('20161130')) +# print('20161201', MgSB.getday('20161201')) +# print('20161202', MgSB.getday('20161202')) +# t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S') +# print(t, MgSB.gettime(t)) plt.figure() MgSB._plot_all_final() - plt.savefig(os.sep.join(('gx','MgSB_plot_all_final.png')), - pad_inches=0.33,bbox_inches='tight',dpi=100) +# plt.savefig(os.sep.join(('gx','MgSB_plot_all_final.png')), +# pad_inches=0.33,bbox_inches='tight',dpi=100) plt.show() diff --git a/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py b/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py index b0c6667..9513114 100644 --- a/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py +++ b/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py @@ -30,7 +30,7 @@ def __init__(self, verbose=False): self._data = {} self.nfinal = 0 first = True - for yyyymmdd in (zM.viewkeys() & zT.viewkeys()): + for yyyymmdd in (zM.keys() & zT.keys()): # load data dictionary rc1, Mg, SB = zM.getday(yyyymmdd) @@ -141,18 +141,18 @@ def output_final_textfile(self, filename): NB = os.environ['NOBACKUP'] z = TSI_Mg_SB_Merged_Daily() - print('16000101', z.getday('16000101')) - print('20161130', z.getday('20161130')) - print('20161201', z.getday('20161201')) - print('20161202', z.getday('20161202')) - t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S') - print(t, z.gettime(t)) +# print('16000101', z.getday('16000101')) +# print('20161130', z.getday('20161130')) +# print('20161201', z.getday('20161201')) +# print('20161202', z.getday('20161202')) +# t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S') +# print(t, z.gettime(t)) z.output_final_textfile(os.sep.join((NB,'NRLSSI2','output','NRLSSI2.txt'))) plt.figure(figsize=(6,12)) z._plot_all_final() - plt.savefig(os.sep.join(('gx','TSInMgSB_plot_all_final.png')), - pad_inches=0.33,bbox_inches='tight',dpi=100) +# plt.savefig(os.sep.join(('gx','TSInMgSB_plot_all_final.png')), +# pad_inches=0.33,bbox_inches='tight',dpi=100) plt.show() diff --git a/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py b/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py index 2f47947..f511378 100644 --- a/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py +++ b/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py @@ -171,8 +171,8 @@ def __init__(self, sys.stdout.write('\n') sys.stdout.flush() - def viewkeys(self): - return self._data.viewkeys() + def keys(self): + return self._data.keys() def getday(self, yyyymmdd): ''' @@ -234,16 +234,16 @@ def _plot_all_final(self): if __name__ == '__main__': TSI = TSI_Daily() - print('16000101', TSI.getday('16000101')) - print('20161130', TSI.getday('20161130')) - print('20161201', TSI.getday('20161201')) - print('20161202', TSI.getday('20161202')) - t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S') - print(t, TSI.gettime(t)) +# print('16000101', TSI.getday('16000101')) +# print('20161130', TSI.getday('20161130')) +# print('20161201', TSI.getday('20161201')) +# print('20161202', TSI.getday('20161202')) +# t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S') +# print(t, TSI.gettime(t)) plt.figure() TSI._plot_all_final() - plt.savefig(os.sep.join(('gx','TSI_plot_all_final.png')), - pad_inches=0.33,bbox_inches='tight',dpi=100) +# plt.savefig(os.sep.join(('gx','TSI_plot_all_final.png')), +# pad_inches=0.33,bbox_inches='tight',dpi=100) plt.show() From 7465c97f149d2a62445659e3b4fafeb6b4e933b0 Mon Sep 17 00:00:00 2001 From: Peter Norris Date: Thu, 4 Apr 2024 10:26:50 -0400 Subject: [PATCH 3/4] pmn: updated and tested code and instructions --- .../NRLSSI2/Mg_SB_from_daily_file.py | 20 +++--- GEOS_RadiationShared/NRLSSI2/README.md | 66 +++++++++++++++++-- .../NRLSSI2/TSI_Mg_SB_merged_from_daily.py | 13 ++-- .../NRLSSI2/TSI_from_daily_files.py | 14 ++-- 4 files changed, 83 insertions(+), 30 deletions(-) diff --git a/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py b/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py index e3735a3..ba9d3a0 100644 --- a/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py +++ b/GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py @@ -20,23 +20,21 @@ def _yyyymmdd(t): class Mg_SB_Daily: ''' Daily time series for Solar Mg and SB indices from NRLSSI2 ''' - def __init__(self, + def __init__(self, DATADIR, filenames=[ - 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20170630_c20170928.txt', - 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20180331_c20180418.txt', - 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20191231_c20200225.txt', - 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20201231_c20210204.txt', - 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20220630_c20220803.txt', - 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20230331_c20230411.txt', - 'data/tsi-ssi_v02r01_model-input-time-series_s18820101_e20231231_c20240221.txt', + 'tsi-ssi_v02r01_model-input-time-series_s18820101_e20170630_c20170928.txt', + 'tsi-ssi_v02r01_model-input-time-series_s18820101_e20180331_c20180418.txt', + 'tsi-ssi_v02r01_model-input-time-series_s18820101_e20191231_c20200225.txt', + 'tsi-ssi_v02r01_model-input-time-series_s18820101_e20201231_c20210204.txt', + 'tsi-ssi_v02r01_model-input-time-series_s18820101_e20220630_c20220803.txt', + 'tsi-ssi_v02r01_model-input-time-series_s18820101_e20230331_c20230411.txt', + 'tsi-ssi_v02r01_model-input-time-series_s18820101_e20231231_c20240221.txt', ], verbose=False): self.filenames = filenames self.verbose = verbose - NB = os.environ['NOBACKUP'] - if verbose: print('reading Mg and SB indicies from', filenames) @@ -45,7 +43,7 @@ def __init__(self, self.nfinal = 0 first = True for filename in filenames: - with open(os.sep.join((NB,'NRLSSI2',filename))) as f: + with open(os.sep.join((DATADIR,filename))) as f: lines = f.readlines() # process the data, ignoring the header (top line) diff --git a/GEOS_RadiationShared/NRLSSI2/README.md b/GEOS_RadiationShared/NRLSSI2/README.md index bc012c2..355e3e0 100644 --- a/GEOS_RadiationShared/NRLSSI2/README.md +++ b/GEOS_RadiationShared/NRLSSI2/README.md @@ -1,14 +1,70 @@ -get data files from https://www.ncei.noaa.gov/data/total-solar-irradiance/access/ -Note: can use wget https://... from discover +========================================================= +Instructions on how to pre-process NRLSSI2 data for GEOS5 +========================================================= +Peter Norris, April 2024 + +A. Introduction. +In the GEOS5 RRTMG and RRTMGP solar codes, the Solar Spectral Irradiance (SSI) at TOA is calculated from +NRLSSI2 Total Solar Irradiance (TSI) and Mg (facular) and SB (sunspot) indices. This data is currently +updated yearly, when the finalized yearly products become available, which is generally in late January +to late February of the following year. For example, the finalized 2023 data came out on 2024-02-23. This +README describes how to get the NRLSSI2 source data and how to pre-process it into a text file read by +the SolarGC. + +B. Get the data for the NEW year. +The source data is obtained from https://www.ncei.noaa.gov/data/total-solar-irradiance/access/ +Note: you can use wget https://... from discover. +There are two data files required for each new year: 1. from ancilliary-data/, e.g., tsi-ssi_v02r01_model-input-time-series_s18820101_e20180331_c20180418.txt !!! NOTE: Please read comments in docstring of Mg_SB_from_daily_file.py BEFORE downloading !!! In particular, - -- the filenames list in that file should ONLY be APPENDED to, + -- the list in that python code should ONLY be APPENDED to, which will ensure that any existing final data is not overwritten, - thus ensuring historical reproduceability. + thus ensuring historical reproducibility. 2. from daily/, e.g., tsi_v02r01_daily_s20170101_e20171231_c20180227.nc !!! NOTE: Please read comments in docstring of TSI_from_daily_files.py BEFORE downloading !!! In particular, -- there should be no files with overlapping time periods, and -- you should NEVER overwrite any existing non-preliminary file, - since this may cause historical non-reproduceability. + since this may cause historical non-reproducibility. +NOTE: Only two files are needed for the new year. Existing data has already been obtained +for previous years. There is no need to get historical data for past years. Even though the +file in part 1 above contains historical data, only the data for the new year inside it will +be used if you follow these instructions, in order to maintain historical reproducibility. + +C. Decide where to put the source data. +The source data is put in the directory DATADIR, where DATADIR is currently + /discover/nobackup/projects/gmao/share/gmao_ops/fvInput/g5gcm/solar/NRLSSI2/data +This directory is managed by the GMAO Operations Team, so if you is not they, you will need to download +the data somewhere else first and then work with "Ops" to add it to that directory. DATADIR here is NOT +an environmental variable, just a shorthand. But whatever DATADIR is used must also be explicity set in +the "main" of the preprocessing driver: TSI_Mg_SB_merged_from_daily.py. + +D. Preparing the preprocessor. +1. Add the filename from B.1. to the filenames list in Mg_SB_from_daily_file.py. Specifically APPEND +it to the END of the list in the default argument of the __init__ of class Mg_SB_Daily +in Mg_SB_from_daily_file.py. +2. Set the OUTDIR in the "main" of TSI_Mg_SB_merged_from_daily.py to where you want to store the +pre-processed file for the neaw year. + +E. Running the preprocessor. +python TSI_Mg_SB_merged_from_daily.py +After you exit the program, rename the NRLSSI2.VYYYY.txt in OUTDIR to the correct year. + +F. Storing the Output file +Get "Ops" to put the new NRLSSI2.VYYYY.txt from your OUTDIR into their directory: + /discover/nobackup/projects/gmao/share/gmao_ops/fvInput/g5gcm/solar + +G. Symlinks +There is a symlink in later "Ops solar" directory named NRLSSI2.txt. +It will typically be pointed to the file NRLSSI2.vYYYY.txt for the new year, but that decision must +be approved by Ops and by the Modelling Team and Data Assimilation Team heads. + +H. Using the data. +The GEOS5 run should have the following lines in its ASGCM.rc: +USE_NRLSSI2: .TRUE. +SOLAR_CYCLE_FILE_NAME: ExtData/g5gcm/solar/NRLSSI2.txt +(Note: These will only be used if the solar code being run is RRTMG or RRTMGP. Chou-Suarez runs do +not use NRLSSI2 data). +PS: You are free to change the use of the symlink NRLSSI2.txt to any specific NRLSSI2.vYYYY.txt +if you want reproducibility to a specific historical run. diff --git a/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py b/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py index 9513114..04aa28d 100644 --- a/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py +++ b/GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py @@ -17,13 +17,13 @@ def _yyyymmdd(t): class TSI_Mg_SB_Merged_Daily: ''' Daily time series for TSI and Solar Mg and SB indices from NRLSSI2 ''' - def __init__(self, verbose=False): + def __init__(self, DATADIR, verbose=False): self.verbose = verbose # get each data set - zM = Mg_SB_Daily(verbose=verbose) - zT = TSI_Daily(verbose=verbose) + zM = Mg_SB_Daily(DATADIR, verbose=verbose) + zT = TSI_Daily(DATADIR, verbose=verbose) # form the INTERSECTION: # (each set is indexed by yyyymmdd, unique in each set) @@ -138,9 +138,10 @@ def output_final_textfile(self, filename): if __name__ == '__main__': - NB = os.environ['NOBACKUP'] + DATADIR = '/discover/nobackup/projects/gmao/share/gmao_ops/fvInput/g5gcm/solar/NRLSSI2/data' + OUTDIR = os.sep.join((os.environ['NOBACKUP'],'NRLSSI2','output')) - z = TSI_Mg_SB_Merged_Daily() + z = TSI_Mg_SB_Merged_Daily(DATADIR) # print('16000101', z.getday('16000101')) # print('20161130', z.getday('20161130')) # print('20161201', z.getday('20161201')) @@ -148,7 +149,7 @@ def output_final_textfile(self, filename): # t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S') # print(t, z.gettime(t)) - z.output_final_textfile(os.sep.join((NB,'NRLSSI2','output','NRLSSI2.txt'))) + z.output_final_textfile(os.sep.join((OUTDIR,'NRLSSI2.vYYYY.txt'))) plt.figure(figsize=(6,12)) z._plot_all_final() diff --git a/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py b/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py index f511378..c83c755 100644 --- a/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py +++ b/GEOS_RadiationShared/NRLSSI2/TSI_from_daily_files.py @@ -1,13 +1,13 @@ ''' Get TSI from daily NRLSSI2 netCDF files ============================================================ -IMPORTANT RULES THAT GOVERN data/tsi_*_daily_*.nc FILES: +IMPORTANT RULES THAT GOVERN tsi_*_daily_*.nc FILES: ============================================================ In order to maintain historical reproducabilty and not fail the error checking below, please ensure the following for -data/tsi_*_daily_*.nc files: +tsi_*_daily_*.nc files: 1. files specified in the file_template argument must -obey 'data/tsi_*_daily_s????????_e????????_c????????.nc' +obey 'tsi_*_daily_s????????_e????????_c????????.nc' 2. the first '*' is the version e.g., 'v02r01', with an optional trailing '-preliminary' 3. The starting date s???????? should be of the form @@ -55,8 +55,8 @@ def _yyyymmdd(t): class TSI_Daily: ''' Daily time series for NRLSSI2 TSI ''' - def __init__(self, - file_template='data/tsi_*_daily_s????????_e????????_c????????.nc', + def __init__(self, DATADIR, + file_template='tsi_*_daily_s????????_e????????_c????????.nc', verbose=False): self.file_template = file_template @@ -69,11 +69,9 @@ def __init__(self, # Validate files # ============== - NB = os.environ['NOBACKUP'] - # collect templated files and validate files = [] - for path in glob(os.sep.join((NB,'NRLSSI2',file_template))): + for path in glob(os.sep.join((DATADIR,file_template))): filename = os.path.basename(path) # validate filename v = filename.split('_') From 158509e17eb35927790b7e9422f47f1a949cdcef Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Mon, 12 Aug 2024 15:10:19 -0400 Subject: [PATCH 4/4] Update ESMF CMake target to ESMF::ESMF This PR updates the ESMF CMake target to `ESMF::ESMF` which is the correct canonical target name for ESMF. This is necessary for Spack compatibility. NOTE: This requires ESMF 8.6.1 or later. --- CMakeLists.txt | 2 +- GEOSsolar_GridComp/RRTMG/rrtmg_sw/gcm_model/src/CMakeLists.txt | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 608f2c4..c1af80e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,5 +10,5 @@ set (alldirs esma_add_library (${this} SRCS GEOS_RadiationGridComp.F90 SUBCOMPONENTS ${alldirs} - DEPENDENCIES MAPL GEOS_Shared esmf) + DEPENDENCIES MAPL GEOS_Shared ESMF::ESMF) diff --git a/GEOSsolar_GridComp/RRTMG/rrtmg_sw/gcm_model/src/CMakeLists.txt b/GEOSsolar_GridComp/RRTMG/rrtmg_sw/gcm_model/src/CMakeLists.txt index 8258262..c93507f 100644 --- a/GEOSsolar_GridComp/RRTMG/rrtmg_sw/gcm_model/src/CMakeLists.txt +++ b/GEOSsolar_GridComp/RRTMG/rrtmg_sw/gcm_model/src/CMakeLists.txt @@ -40,7 +40,7 @@ foreach (file ${k_g_srcs}) endforeach () esma_add_library (${this} SRCS ${srcs} ${k_g_srcs} - DEPENDENCIES GEOS_Shared MAPL RRTMG_SW_mods GEOS_RadiationShared esmf NetCDF::NetCDF_Fortran) + DEPENDENCIES GEOS_Shared MAPL RRTMG_SW_mods GEOS_RadiationShared ESMF::ESMF NetCDF::NetCDF_Fortran) if(ENABLE_SOLAR_RADVAL) target_compile_definitions(${this} PRIVATE SOLAR_RADVAL)