In [1]:
%run basics
import pytz
import re
%matplotlib
In [2]:
def get_clearsky_Fsd(ddd,hhh,k,info):
longitude = info["longitude"]
latitude = info["latitude"]
altitude = info["altitude"]
GMT_zone = info["GMT_zone"]
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-longitude)/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(latitude))*numpy.sin(declination)+
numpy.cos(numpy.radians(latitude))*numpy.cos(declination)*
numpy.cos(hour_angle))
zenith = numpy.ma.masked_greater_equal(zenith,numpy.pi/2)
air_mass_ratio = (numpy.exp(-1*altitude/8343.5)/
(numpy.cos(zenith)+0.15*(numpy.degrees(90-zenith)+3.855)**-1.253))
Fsd_cs = Fsd_TOA*numpy.exp(-k*air_mass_ratio)*numpy.cos(zenith)
return Fsd_cs
In [3]:
def residuals(params,Fsd_obs,ddd,hhh,info,weighted,weight_type):
k = params[0]
nDays = info["nDays"]
ntsInDay = info["ntsInDay"]
Fsd_obs_max = numpy.ma.max(Fsd_obs.reshape(nDays,ntsInDay),axis=0)
Fsd_cs = get_clearsky_Fsd(ddd,hhh,k,info)
Fsd_cs_max = numpy.ma.max(Fsd_cs.reshape(nDays,ntsInDay),axis=0)
resid = Fsd_obs_max - Fsd_cs_max
if weighted:
weights = CalculateWeights(resid,weight_type=weight_type)
return resid*weights
else:
return resid
In [4]:
def irls_leastsq(func,p0,args=(),maxfev=3,weight_type='Huber',mode='verbose'):
weighted = False
popt,pcov = scipy.optimize.leastsq(func,p0,args=args+(weighted,weight_type))
if mode!="quiet":
print 'After non-weighted call to leastsq: ',popt
n=1
weighted = True
while n<=maxfev:
popt,pcov = scipy.optimize.leastsq(func,popt,args=args+(weighted,weight_type))
if mode!="quiet":
print 'Weighted call '+str(n)+' to leastsq: ',popt
n = n + 1
return popt,pcov
def CalculateWeights(resid,weight_type='Huber'):
if weight_type.lower()=='tukey':
weights = TukeyBiweight(MAD_scale(resid))
elif weight_type.lower()=='huber':
weights = HuberTWeight(MAD_scale(resid))
elif weight_type.lower()=='ols':
weights = OLSWeight(MAD_scale(resid))
else:
print "CalculateWeights: Unknown weight type, used Huber's t weight"
weights = HuberTWeight(MAD_scale(resid))
return weights
def TukeyBiweight(x):
return ((numpy.abs(x)<1.0)*(1.0 - x**2)**2)
def HuberTWeight(x):
weights = numpy.ones(len(x))
idx = numpy.where(numpy.abs(x)>1.345)[0]
weights[idx] = 1.345/numpy.abs(x[idx])
return weights
def OLSWeight(x):
weights = numpy.ones(len(x))
return weights
def MAD_scale(resid):
scale = numpy.median(numpy.abs(resid - numpy.median(resid))) / 0.6745
return resid/scale
In [19]:
# read the netCDF file
ncname = qcio.get_filename_dialog()
ds = qcio.nc_read_series(ncname)
dt = ds.series["DateTime"]["Data"]
In [20]:
# create the info dictionary to hold useful stuff
info = {}
ts = int(ds.globalattributes["time_step"])
info["ntsInDay"] = int(24.0*60.0/float(ts))
info["longitude"] = float(ds.globalattributes["longitude"])
info["latitude"] = float(ds.globalattributes["latitude"])
if "elevation" in ds.globalattributes:
info["altitude"] = int(re.sub("[^0123456789\.]","",str(ds.globalattributes["elevation"])))
else:
info["altitude"] = 1.0
time_zone = ds.globalattributes["time_zone"]
tz = pytz.timezone(time_zone)
cdt = datetime.datetime(2013,6,1,0,0,0)
utcoffset = tz.utcoffset(cdt)
info["GMT_zone"] = utcoffset.seconds/float(3600)
In [21]:
# do each year and month separately
start_date = dt[0]
end_date = start_date+dateutil.relativedelta.relativedelta(months=1)
last_date = dt[-1]
while start_date<last_date:
si = qcutils.GetDateIndex(dt,str(start_date),ts=ts,match="startnextday",default=0)
ei = qcutils.GetDateIndex(dt,str(end_date),ts=ts,match="endpreviousday",default=len(dt)-1)
mi = int((si+ei)/2)
info["nDays"] = int((ei-si+1)/info["ntsInDay"])
ddd,f,a = qcutils.GetSeries(ds,"Ddd",si=si,ei=ei)
hhh,f,a = qcutils.GetSeries(ds,"Hdh",si=si,ei=ei)
Fsd,f,a = qcutils.GetSeries(ds,"Fsd",si=si,ei=ei)
# correct decimal hour to refer to centre of period rather than end of period
hhh = hhh - float(ts)/(60.0*2)
popt,pcov = irls_leastsq(residuals,[0.27],args=(Fsd,ddd,hhh,info),maxfev=3,weight_type='huber',mode="quiet")
print dt[si],dt[mi],dt[ei],popt[0]
start_date = end_date
end_date = start_date+dateutil.relativedelta.relativedelta(months=1)
In [22]:
month,f,a = qcutils.GetSeries(ds,"Month")
ddd_all,f,a = qcutils.GetSeries(ds,"Ddd")
hhh_all,f,a = qcutils.GetSeries(ds,"Hdh")
Fsd_all,f,a = qcutils.GetSeries(ds,"Fsd")
for m in range(1,13):
idx = numpy.where(month==m)[0]
ddd = ddd_all[idx]
hhh = hhh_all[idx]
Fsd = Fsd_all[idx]
info["nDays"] = int(len(idx)/info["ntsInDay"])
# correct decimal hour to refer to centre of period rather than end of period
hhh = hhh - float(ts)/(60.0*2)
popt,pcov = irls_leastsq(residuals,[0.27],args=(Fsd,ddd,hhh,info),maxfev=3,weight_type='huber',mode="quiet")
print m,popt[0]
In [ ]: