In [1]:
%run basics
import pytz
import re
%matplotlib
In [2]:
def get_clearsky_Fsd(ds):
"""
Calculate downwelling radiation at surface.
Code from Ian McHugh.
"""
# set a value for "k", the extinction coefficient
k = 0.25
# get the latitude and longitude
lat_decdeg = float(ds.globalattributes["latitude"])
long_decdeg = float(ds.globalattributes["longitude"])
if "elevation" in ds.globalattributes:
ALT_m = int(re.sub("[^0123456789\.]","",str(ds.globalattributes["elevation"])))
else:
ALT_m = 1.0
rec_length = int(ds.globalattributes["time_step"])
# Ian's original code queried AskGeo for the UTC offset in milliseconds
# and converted that to hours. Here we use the time zone given in the
# netCDF file global attributes.
# get the time zone
time_zone = ds.globalattributes["time_zone"]
# get a pytz object
tz = pytz.timezone(time_zone)
# use the same "current" date as Ian
cdt = datetime.datetime(2013,6,1,0,0,0)
# get the offset to UTC time as a timedelta object
utcoffset = tz.utcoffset(cdt)
# get the offset as hours
GMT_zone = utcoffset.seconds/float(3600)
# get the datetime
ldt = ds.series["DateTime"]["Data"]
# get the day of the year from local time
#DOY = numpy.array([dt.timetuple().tm_yday for dt in ldt])
start_doy = ldt[0].timetuple().tm_yday
end_doy = ldt[-2].timetuple().tm_yday
DOY = numpy.arange(start_doy,end_doy+1,1)
# For each day calculate equation of time correction, solar noon, declination and TOA radiation
array_EqofTime=0.17*numpy.sin(4*numpy.pi*(DOY-80)/373)-0.129*numpy.sin(2*numpy.pi*(DOY-8)/355) # DiLaura (1984)
array_solar_noon=12+(GMT_zone*15.0-long_decdeg)/360*24-array_EqofTime # Me
array_decl=numpy.radians(23.4)*numpy.sin((DOY+284)/365.0*2*numpy.pi) # Oke (1987)
array_TOArad=(1+0.034*numpy.cos(DOY/365.25*2*numpy.pi))*1367.0 # Duffie and Beckman (1980)
# Create an hour angle array for each minute of day and each day of year
array_h=numpy.tile(numpy.linspace(0,1439.0/1440*24,num=1440),(len(DOY),1))
array_h=abs(numpy.radians((array_solar_noon.reshape(len(DOY),1)-array_h)*15))
# Duplicate declination array for each time of day
array_decl=numpy.tile(array_decl,(1440,1)).T
# Calculate zenith angles
array_z=numpy.arccos(numpy.sin(numpy.radians(lat_decdeg))*numpy.sin(array_decl)+numpy.cos(numpy.radians(lat_decdeg))*numpy.cos(array_decl)*numpy.cos(array_h))
array_z_msk=numpy.ma.masked_greater_equal(array_z,numpy.pi/2) # Mask night values
# Calculate optical air mass term for all valid Z
array_m=(numpy.exp(-1*ALT_m/8343.5)/(numpy.cos(array_z_msk)+0.15*(numpy.degrees(90-array_z_msk)+3.855)**-1.253)) # Wunderlich (1972)
# Instantaneous clear sky surface radiation in Wm-2 for each minute of the day
array_Kdown_clr_mins=array_TOArad.reshape(len(array_TOArad),1)*numpy.exp(-k*array_m)*numpy.cos(array_z_msk)
# Aggregate one-minute instantaneous clear sky rad to period average
array_Kdown_clr_hr=numpy.empty([len(DOY),1440/rec_length])
for i in xrange(len(DOY)):
array_temp=array_Kdown_clr_mins[i][:].reshape(1440/rec_length,rec_length) # Temporary bins
array_Kdown_clr_hr[i][:]=numpy.ma.mean(array_temp,axis=1) # Average bin content
# Aggregate to daily
array_Kdown_clr_daily=(array_Kdown_clr_hr*(rec_length*60.0/10**6)).sum(axis=1)
#return array_Kdown_clr_daily,array_Kdown_clr_hr
return array_Kdown_clr_hr.flatten(order="C")
In [ ]:
def get_clearsky_radiation(doy,params):
longitude = params["longitude"]
GMT_zone = params["GMT_zone"]
eot = equation_of_time(doy)
solar_noon = get_solar_noon(eot,longitude,GMT_zone)
declination = get_declination(doy)
flux_toa = get_radiation_topofatmosphere(doy)
In [3]:
# read the netCDF file
ncname = qcio.get_filename_dialog()
ds = qcio.nc_read_series(ncname)
ts = int(ds.globalattributes["time_step"])
ntsInDay = int(24.0*60.0/float(ts))
In [4]:
# get the clear sky radiation
Fsd_cs = get_clearsky_Fsd(ds)
In [5]:
# get the measured downwelling shortwave
Fsd,f,a = qcutils.GetSeriesasMA(ds,"Fsd")
In [6]:
# plot
ldt = ds.series["DateTime"]["Data"]
fig = plt.figure()
plt.plot(ldt,Fsd,'b.')
plt.plot(ldt,Fsd_cs,'r+')
plt.show()
In [7]:
start_date = ldt[0]
end_date = start_date+dateutil.relativedelta.relativedelta(months=1)
si = qcutils.GetDateIndex(ldt,str(start_date),ts=ts,match="startnextday",default=0)
ei = qcutils.GetDateIndex(ldt,str(end_date),ts=ts,match="endpreviousday",default=len(ldt))
nDays = int((ei-si+1)/ntsInDay)
Fsd_period = Fsd[si:ei+1].reshape(nDays,ntsInDay)
Fsd_cs_period = Fsd_cs[si:ei+1].reshape(nDays,ntsInDay)
Fsd_max = numpy.ma.max(Fsd_period,axis=0)
Fsd_cs_max = numpy.ma.max(Fsd_cs_period,axis=0)
In [8]:
fig = plt.figure()
plt.plot(numpy.arange(0,len(Fsd_max)),Fsd_max,'b.')
plt.plot(numpy.arange(0,len(Fsd_max)),Fsd_cs_max,'r+')
plt.show()
In [29]:
k = 0.25
time_zone = ds.globalattributes["time_zone"]
long_decdeg = float(ds.globalattributes["longitude"])
lat_decdeg = float(ds.globalattributes["latitude"])
if "elevation" in ds.globalattributes:
ALT_m = int(re.sub("[^0123456789\.]","",str(ds.globalattributes["elevation"])))
else:
ALT_m = 1.0
rec_length = int(ds.globalattributes["time_step"])
tz = pytz.timezone(time_zone)
cdt = datetime.datetime(2013,6,1,0,0,0)
utcoffset = tz.utcoffset(cdt)
GMT_zone = utcoffset.seconds/float(3600)
start_doy = ldt[0].timetuple().tm_yday
end_doy = ldt[-2].timetuple().tm_yday
DOY = numpy.arange(start_doy,end_doy+1,1)
array_EqofTime=0.17*numpy.sin(4*numpy.pi*(DOY-80)/373)-0.129*numpy.sin(2*numpy.pi*(DOY-8)/355)
array_solar_noon=12+(GMT_zone*15.0-long_decdeg)/360*24-array_EqofTime
array_decl=numpy.radians(23.4)*numpy.sin((DOY+284)/365.0*2*numpy.pi)
array_TOArad=(1+0.034*numpy.cos(DOY/365.25*2*numpy.pi))*1367.0
array_hh=numpy.tile(numpy.linspace(0,1439.0/1440*24,num=1440),(len(DOY),1))
array_h=abs(numpy.radians((array_solar_noon.reshape(len(DOY),1)-array_hh)*15))
array_decl=numpy.tile(array_decl,(1440,1)).T
array_z=numpy.arccos(numpy.sin(numpy.radians(lat_decdeg))*numpy.sin(array_decl)+numpy.cos(numpy.radians(lat_decdeg))*numpy.cos(array_decl)*numpy.cos(array_h))
array_z_msk=numpy.ma.masked_greater_equal(array_z,numpy.pi/2)
array_m=(numpy.exp(-1*ALT_m/8343.5)/(numpy.cos(array_z_msk)+0.15*(numpy.degrees(90-array_z_msk)+3.855)**-1.253))
array_Kdown_clr_mins=array_TOArad.reshape(len(array_TOArad),1)*numpy.exp(-k*array_m)*numpy.cos(array_z_msk)
array_Kdown_clr_hr=numpy.empty([len(DOY),1440/rec_length])
for i in xrange(len(DOY)):
array_temp=array_Kdown_clr_mins[i][:].reshape(1440/rec_length,rec_length) # Temporary bins
array_Kdown_clr_hr[i][:]=numpy.ma.mean(array_temp,axis=1) # Average bin content
Fsd_cs_im = array_Kdown_clr_hr.flatten(order="C")
In [69]:
k = 0.208
ddd,f,a = qcutils.GetSeries(ds,"Ddd")
hhh,f,a = qcutils.GetSeries(ds,"Hdh")
# correct decimal hour to refer to centre of period rather than end of period
hhh = hhh - float(ts)/(60.0*2)
eot = 0.17*numpy.sin(4*numpy.pi*(ddd-80)/373)-0.129*numpy.sin(2*numpy.pi*(ddd-8)/355)
solar_noon = 12+(GMT_zone*15.0-long_decdeg)/360*24-eot
declination = numpy.radians(23.4)*numpy.sin((ddd+284)/365.0*2*numpy.pi)
Fsd_TOA = (1+0.034*numpy.cos(ddd/365.25*2*numpy.pi))*1367.0
hour_angle = abs(numpy.radians((solar_noon-hhh)*15))
zenith = numpy.arccos(numpy.sin(numpy.radians(lat_decdeg))*numpy.sin(declination)+numpy.cos(numpy.radians(lat_decdeg))*numpy.cos(declination)*numpy.cos(hour_angle))
zenith = numpy.ma.masked_greater_equal(zenith,numpy.pi/2)
air_mass_ratio = (numpy.exp(-1*ALT_m/8343.5)/(numpy.cos(zenith)+0.15*(numpy.degrees(90-zenith)+3.855)**-1.253))
Fsd_cs_pri = Fsd_TOA*numpy.exp(-k*air_mass_ratio)*numpy.cos(zenith)
In [17]:
fig = plt.figure()
plt.plot(array_hh[0,:],array_h[0,:])
plt.plot(hhh[0:47],hour_angle[0:47],'r+')
plt.show()
In [24]:
fig2 = plt.figure()
plt.plot(array_hh[0,:],array_z_msk[0,:])
plt.plot(hhh[0:47],zenith[0:47],'r+')
plt.show()
In [28]:
fig3 = plt.figure()
plt.plot(array_hh[0,:],array_m[0,:])
plt.plot(hhh[0:47],air_mass_ratio[0:47],'r+')
plt.show()
In [36]:
fig4 = plt.figure()
plt.plot(hhh[0:47],Fsd_cs_im[0:47])
plt.plot(hhh[0:47],Fsd_cs_pri[0:47],'r+')
plt.show()
In [70]:
Fsd_cs_im_period = Fsd_cs_im[si:ei+1].reshape(nDays,ntsInDay)
Fsd_cs_im_max = numpy.ma.max(Fsd_cs_im_period,axis=0)
Fsd_cs_pri_period = Fsd_cs_pri[si:ei+1].reshape(nDays,ntsInDay)
Fsd_cs_pri_max = numpy.ma.max(Fsd_cs_pri_period,axis=0)
In [71]:
fig5 = plt.figure()
plt.plot(numpy.arange(0,len(Fsd_max)),Fsd_max,'b.')
plt.plot(numpy.arange(0,len(Fsd_max)),Fsd_cs_im_max,'r+')
plt.plot(numpy.arange(0,len(Fsd_max)),Fsd_cs_pri_max,'g^')
plt.show()
In [ ]: