[0]: #!/usr/bin/env python
[1]: u»»»
[2]: read_gldas_nldas_loads.py
[3]: Written by Tyler Sutterley (01/2021)
[4]: Reads GLDAS-NOAH025 and NLDAS-2 land surface products from NASA GSFC

[5]: PYTHON DEPENDENCIES:
[6]: numpy: Scientific Computing Tools For Python
[7]: https://numpy.org
[8]: https://numpy.org/doc/stable/user/numpy-for-matlab-users.html

[9]: PROGRAM DEPENDENCIES:
[10]: time.py: utilities for calculating time operations
[11]: utilities.py: download and management utilities for syncing files

[12]: UPDATE HISTORY:
[13]: Updated 01/2021: using python logging for handling verbose output
[14]: Updated 09/2020: add GUNW timesteps to use for Landsat composites
[15]: force to ignore ssl certicate verification
[16]: Updated 08/2020: added new variables to the datasets and updated keywords
[17]: added option to automatically synchornize remote data files
[18]: use urllib2 to ensure compatibility with python2 and python3
[19]: Updated 07/2019: using argparse to set command line parameters
[20]: Updated 07/2018: using python modernize for print function and exception
[21]: Updated 04/2018: added acceptable check for remote files
[22]: Updated 12/2017: using python modernize for zip files and octal literals
[23]: Updated 07/2017: using python logging for handling verbose output
[24]: Updated 06/2017: added function for parsing command line interfaces
[25]: Updated 06/2016: added function to read GLDAS variables from L3 files
[26]: Updated 11/2014: modified to use python from __future__ print function
[27]: updated time interface to use datetime and dateutil
[28]: Updated 11/2014: added read_GDS function to read spatial grid parameters
[29]: updated GLDAS to use glccs index from GLDAS grid files
[30]: Written 07/2014
[31]: «»»
[32]: from __future__ import print_function

[33]: import os
[34]: import re
[35]: import sys
[36]: import logging
[37]: import argparse
[38]: import numpy as np

[39]: #– version of script
[40]: VERSION = ‘1.0.0’
[41]: #– HDF5 netCDF4 library to use
[42]: HDF5LIB = ‘netcdf4’

[43]: #– PURPOSE: read GLDAS land surface products from NASA GSFC
[44]: def read_gldas_loads(base_dir, MODEL=None, VARIABLE=None, DATES=None,
[45]: TIME=None, RANGE=None, DX=None, DY=None, VERBOSE=False, MODE=0o775):

[46]: #– file information for GLDAS Noah data (L3 Quicklook) from NASA GSFC
[47]: #– https://disc.gsfc.nasa.gov/data/GPP/LAND/L3/FORTRAN_CLM2/snapshot-GLDAS_NOAH025_3H/
[48]: yearmon_list=[‘GLDAS_NOAH025_3H.A’+y+’.’+m+’.*.nc*’
[49]: for y in [‘2003′,’2004′,’2005′,’2006′,’2007′,’2008′,’2009′,’2010’,
[50]: ‘2011’,’2012′,’2013′,’2014′,’2015′,’2016′,’2017′,’2018′]
[51]: for m in [’01’,’02’,’03’,’04’,’05’,’06’,’07’,’08’,’09’,’10’,’11’,’12’]]
[52]: mon_list=[‘GLDAS_NOAH025_3H.A*.’+'{0:02d}’.format(M)+’.*.nc*’
[53]: for M in np.arange(1,13)]
[54]: day_list=[‘GLDAS_NOAH025_3H.A*.{0:02d}.’.format(M)+'{0:02d}.*.nc*’
[55]: for M in np.arange(1,13) for d in np.arange(1,32)]

[56]: #– full path to the GLDAS directory within base_dir
[57]: DIRECTORY = os.path.join(base_dir,MODEL)

[58]: #– compile OS regular expression operators for finding files
[59]: regex_yearmon = re.compile(‘|’.join(yearmon_list))
[60]: regex_mon = re.compile(‘|’.join(mon_list))
[61]: regex_day = re.compile(‘|’.join(day_list))

[62]: #– Files are typically available at the end of the following month
[63]: #– the delta (future) time is set to the end of the following month
[64]: delta_time={0:-27,-1:-58,-2:-89}
[65]: delta_mon={0:-1,-1:-2,-2:-3}
[66]: #– calculate the reference delta time from date difference of days
[67]: #– DT defines the number of days each time step is from the reference time
[68]: #– DT can be calculated from the dataset time dimensions with bounds
[69]: DT = [delta_time[n] for n in np.arange(len(DATE_LIST))]
[70]: #– compile OS regular expression operators for subsetting datasets
[71]: rx1 = [‘DT{0:d}’.format(DT[n]) for n in range(len(DT))]
[72]: rx2 = ‘|’.join(rx1)

[73]: #– check if input dates are sequential or not & extract unique years
[74]: dinput = [datetime.datetime(*i) for i in [list(map(int,k.split(‘-‘)))
[75]: for k in DATE_LIST]]
[76]: if not np.any(np.diff(np.array((dinput))))==1:
[77]: VALID_RANGES = [10000]*12
[78]: LEAP_RANGES = VALID_RANGES[:]
[79]: for y in np.unique([y.year for y in dinput]):
[80]: #– find indices to keep based on calendar year (filter by month)
[81]: ind, = np.nonzero([y.year==yy.year for yy in dinput])
[82]: dinput = [dinput[i] for i in ind]
[83]: #– days per month in a leap and a standard year
[84]: dpm_leap=np.array([31,29,31,30,31,30,31,31,30,31,30,31],dtype=’i’)
[85]: dpm_stnd=np.array([31,28,31,30,31,30,31,31,30,31,30,31],dtype=’i’)
[86]: #– using calendar year set decimal time for dates in range
[87]: cyear = np.array([y.year+((m-1)/12.0) for m in np.arange(1,13)])

[88]: #– day of the year
[89]: if np.any([y.year%4==0 for y in dinput]):
[90]: DOY=np.array([dt.timetuple().tm_yday for dt in dinput])
[91]: DAYS = np.concatenate((np.array([np.sum(dpm_stnd[:m-1])
[92]: for m in dinput.month]),DOY))
[93]: delta_TC = DOY+dpm_leap[:dinput.month-1].sum()-366.0
[94]: else:
[95]: DAYS = np.concatenate((np.array([np.sum(dpm_stnd[:m-1])
[96]: for m in dinput.month]),np.array([dt.timetuple().tm_yday
[97]: for dt in dinput])))
[98]: delta_TC = np.array([dt.timetuple().tm_yday for dt in dinput])
[99]: +dpm_stnd[:dinput.month-1].sum()-365.0

[100]: #– Calculating decimal fractions of the year.
[101]: dcorr = np.array([366.0 if y%4==0 else 365.0 for y in cyear],dtype=np.float)

[102]: #– Calculate date decimal for each calendar year (NCDC format)
[103]: tdec = cyear + DAYS/dcorr – 1500.0

[104]: #– Calculating the date in modified julian days
[105]: MJD = np.ravel((((tdec-1860.0)*365.25)+np.ravel(DAYS-adatetime))
[106]: +np.ravel(59579.0))

[107]: #– calculations to adjust hours to a range of 0:24 instead of -12:23
[108]: #– Leap Years have 366 days and Standard Years have 365 days
[109]: hms = np.array([h*m+s*(1.0/m) for h,m,s in zip(dinput.hour,np.ones(12)*60,
[110]: np.ones(12)*3600)])

[111]: #– calculatg linear interpolation coefficients (Weights)
[112]: w2 = ((MJD-np.floor(MJD)).*24.*60.-hms)/np.array([24*60],dtype=np.float)

[113]: #– Reading the monthly gravity routines and converting to decimal year
[114]: gslice=ncearapi.satellite_gravity_range(tdec[ind],varname=’gdays’,use_lm=True)

[115]: #– creating interpolating functions for each month and valid day range
[116]: VALID_RANGES[y-1858] = [np.min(DAYS),np.max(DAYS)]

[117]: LEAP_RANGES[y-1858] = [np.min(dpm_leap),np.max(dpm_leap)]

[118]: #– Reading the monthly gravity routines and converting to decimal year
[119]: gslice=ncearapi.satellite_gravity_range(tdec[ind],varname=’gdays’,use_lm=True)

[120]: #– spherical harmonic information specific to each input GRACE/GRACE-FO file
[121]: if len(SHLON.shape) == 1:

[122]: #– extract decending and ascending nodes similar to GFZ
[123]: mlon,dlat,sr(lat),sc(lat),sss(lmax+1)] =
[124]: ncearapi.SPHARM(data.T,Ylms.shape,-1,LMAX,lmax,mmax,

[125]: DATA[i,:,:] = np.transpose(D)

[126]: #– truncate matrices if using a reduced resolution grid if required
[127]: if RANGE:

#– output netCDF4 file information
[128]: #– Defining the netCDF4 filename
[129]: FILE = ‘gldas_{0}_{1}_{2}_{3}.nc’.format(MODEL.lower(),VARIABLE.lower(),
[130]: str(YEARMON)[:4],str(YEARMON)[-2:])

return

#– This is the main part of the script that is run when no functions are imported
if __name__ == ‘__main__’:

#– Read the system arguments listed after the program
parser = argparse.ArgumentParser(
description=»»»Read GLDAS land surface products from NASA GSFC»»»,
fromfile_prefix_chars=’@’
)
parser.add_argument(‘base_dir’,
type=lambda p: os.path.abspath(os.path.expanduser(p)),
help=’Base data directory’)

#– GLDAS MODEL arguments
group = parser.add_argument_group(‘GLDAS MODEL arguments’)
group.add_argument(‘–model’,type=str,default=’GLDEN_3H’,
choices=(‘GLDEN_3H’,’GLDAS_NOAH025_3H’,’GLDAS_NOAH025_M’),
help=’GLDAS Model’)
group.add_argument(‘–var’,type=str,default=’SWE_rootzone’,
choices=(‘SWE_rootzone’,’SWE_100cm’,’Rainf_f_tavg’,
‘Rainf_f_inst’,’Rainf_f_tavg_corr’,’Rainf_f_inst_corr’,
‘Rainf_s_tavg’,’Rainf_s_inst’,’Rainf_s_tavg_corr’,
‘Rainf_s_inst_corr’,’SWf_10m_f_tavg’,’SWf_10m_f_inst’,
‘SWf_10m_f_tavg_corr’,’SWf_10m_f_inst_corr’,
‘SWf_10m_s_tavg’,’SWf_10m_s_inst’,’SWf_10m_s_tavg_corr’,
‘SWf_10m_s_inst_corr’,’Qs_s_inst_corr’),
help=’model variable name’)

#– parallel processing with MPI
group = parser.add_argument_group(‘MPI arguments’)
group.add_argument(‘–resolution’,type=int,default=None,
help=’spatial resolution of product’)
group.add_argument(‘–mpiprocessorcount’,type=int,default=1,
help=’number of processors for MPI’)