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 [14]:
# 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 [15]:
# 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 [12]:
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 [13]:
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 [17]:
import re
s = "64.5m"
print re.sub("[^0123456789\.]","",s)
In [ ]: