In [1]:
%run basics
%matplotlib
import pysolar
import pytz
from scipy.interpolate import InterpolatedUnivariateSpline
In [2]:
cf = qcio.load_controlfile(path="controlfiles/")
In [ ]:
erai_name = os.path.join(cf["Files"]["erai_path"],cf["Files"]["erai_file"])
erai_file = netCDF4.Dataset(erai_name)
latitude = erai_file.variables["latitude"][:]
longitude = erai_file.variables["longitude"][:]
lat_resolution = abs(latitude[-1]-latitude[0])/(len(latitude)-1)
lon_resolution = abs(longitude[-1]-longitude[0])/(len(longitude)-1)
# get the time and convert to Python datetime object
erai_time = erai_file.variables["time"][:]
time_units = getattr(erai_file.variables["time"],"units")
dt_erai = netCDF4.num2date(erai_time,time_units)
hour_utc = numpy.array([dt.hour for dt in dt_erai])
erai_timestep = 180
In [ ]:
# get a list of sites
site_list = cf["Sites"].keys()
In [ ]:
# now loop over the sies
for site in site_list:
# get the metadata from the control file
out_filename = os.path.join(cf["Sites"][site]["out_filepath"],
cf["Sites"][site]["out_filename"])
site_name = cf["Sites"][site]["site_name"]
site_timezone = cf["Sites"][site]["site_timezone"]
site_latitude = cf["Sites"][site]["site_latitude"]
site_longitude = cf["Sites"][site]["site_longitude"]
site_timestep = cf["Sites"][site]["site_timestep"]
# index of the site in latitude dimension
site_lat_index = int(((latitude[0]-site_latitude)/lat_resolution)+0.5)
erai_latitude = latitude[site_lat_index]
# index of the site in longitude dimension
site_lon_index = int(((site_longitude-longitude[0])/lon_resolution)+0.5)
erai_longitude = longitude[site_lon_index]
# get the UTC and local datetime series
site_tz = pytz.timezone(site_timezone)
# make utc_dt timezone aware
dt_erai_utc = [x.replace(tzinfo=pytz.utc) for x in dt_erai]
# get local time from UTC
dt_erai_loc = [x.astimezone(site_tz) for x in dt_erai_utc]
#dt_erai_loc = [x.astimezone(site_tz) for x in dt_erai_utc]
# NOTE: will have to disable daylight saving at some stage, towers stay on Standard Time
# PRI hopes that the following line will do this ...
#dt_ecmwf_loc = [x-x.dst() for x in dt_ecmwf_loc]
# make local time timezone naive to match datetimes in OzFluxQC
dt_erai_loc_ntz = [x.replace(tzinfo=None) for x in dt_erai_loc]
dt_erai_utc_ntz = [x.replace(tzinfo=None) for x in dt_erai_utc]
# get the datetime in the middle of the accumulation period
erai_offset = datetime.timedelta(minutes=float(erai_timestep)/2)
# subtract this from the ERA-I datetime to align ERA-I 3 hourly points with the middle
# of the accumulation period
dt_erai_loc_cor = [x - erai_offset for x in dt_erai_loc_ntz]
dt_erai_utc_cor = [x - erai_offset for x in dt_erai_utc_ntz]
# now we get the datetime series at the tower time step
tdts = datetime.timedelta(minutes=tower_timestep)
# local time for plotting
start_date = dt_erai_loc_cor[0]
end_date = dt_erai_loc_cor[-1]
dt_erai_loc_tts = [x for x in qcutils.perdelta(start_date,end_date,tdts)]
# UTC for calculating solar altitude
start_date = dt_erai_utc_cor[0]
end_date = dt_erai_utc_cor[-1]
dt_erai_utc_tts = [x for x in qcutils.perdelta(start_date,end_date,tdts)]
# UTC at tower time step as number for interpolation
erai_time_tts = netCDF4.date2num(dt_erai_utc_tts,time_units)
erai_time_3hr = netCDF4.date2num(dt_erai_utc_cor,time_units)
# get the solar altitude, we will use this later to interpolate the ERA Interim
# data from the ERA-I 3 hour time step to the tower time step.
# NOTE: alt_solar is in degrees
alt_solar_3hr = numpy.array([pysolar.GetAltitude(erai_latitude,erai_longitude,dt) for dt in dt_erai_utc_cor])
# get the solar altitude at the tower time step
alt_solar_tts = numpy.array([pysolar.GetAltitude(erai_latitude,erai_longitude,dt) for dt in dt_erai_utc_tts])
idx = numpy.where(alt_solar_tts<=0)[0]
alt_solar_tts[idx] = float(0)
# Interpolate the 3 hourly accumulated downwelling shortwave to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Fsd_3d = erai_file.variables["ssrd"][:,:,:]
Fsd_accum = Fsd_3d[:,site_lat_index,site_lon_index]
# Downwelling shortwave in ERA-I is a cummulative value that is reset to 0 at 0300 and 1500 UTC.
# Here we convert the cummulative values to 3 hourly values.
Fsd_erai_3hr = numpy.ediff1d(Fsd_accum,to_begin=0)
idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]
Fsd_erai_3hr[idx] = Fsd_accum[idx]
Fsd_erai_3hr = Fsd_erai_3hr/(erai_timestep*60)
# normalise the ERA-I downwelling shortwave by the solar altitude
coef_3hr = Fsd_erai_3hr/numpy.sin(numpy.deg2rad(alt_solar_3hr))
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, coef_3hr, k=1)
# get the coefficient at the tower time step
coef_tts = s(erai_time_tts)
# get the downwelling solar radiation at the tower time step
Fsd_erai_tts = coef_tts*numpy.sin(numpy.deg2rad(alt_solar_tts))