In [1]:
%run basics
import pytz
import math
%matplotlib
In [25]:
def get_synthetic_fsd_im(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(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
In [28]:
def GetRadiationDirect(utc_datetime, altitude_deg):
# from Masters, p. 412
if(altitude_deg > 0):
day = GetDayOfYear(utc_datetime)
flux = GetApparentExtraterrestrialFlux(day)
optical_depth = GetOpticalDepth(day)
air_mass_ratio = GetAirMassRatio(altitude_deg)
#return flux * math.exp(-1 * optical_depth * air_mass_ratio)
return flux * math.exp(-0.25 * optical_depth * air_mass_ratio)*math.sin(math.radians(altitude_deg))
else:
return 0.0
In [4]:
def GetAltitudeFast(latitude_deg, longitude_deg, utc_datetime):
# expect 19 degrees for solar.GetAltitude(42.364908,-71.112828,datetime.datetime(2007, 2, 18, 20, 13, 1, 130320))
day = GetDayOfYear(utc_datetime)
declination_rad = math.radians(GetDeclination(day))
latitude_rad = math.radians(latitude_deg)
hour_angle = GetHourAngle(utc_datetime, longitude_deg)
first_term = math.cos(latitude_rad) * math.cos(declination_rad) * math.cos(math.radians(hour_angle))
second_term = math.sin(latitude_rad) * math.sin(declination_rad)
return math.degrees(math.asin(first_term + second_term))
In [5]:
def GetDayOfYear(utc_datetime):
year_start = datetime.datetime(utc_datetime.year, 1, 1, tzinfo=utc_datetime.tzinfo)
delta = (utc_datetime - year_start)
return delta.days
In [6]:
def GetApparentExtraterrestrialFlux(day):
# from Masters, p. 412
return 1160 + (75 * math.sin(math.radians((360./365) * (day - 275))))
In [7]:
def GetOpticalDepth(day):
# from Masters, p. 412
return 0.174 + (0.035 * math.sin(math.radians((360./365) * (day - 100))))
In [8]:
def GetAirMassRatio(altitude_deg):
# from Masters, p. 412
# warning: pukes on input of zero
return (1/math.sin(math.radians(altitude_deg)))
In [9]:
def GetDeclination(day):
'''The declination of the sun is the angle between
Earth's equatorial plane and a line between the Earth and the sun.
The declination of the sun varies between 23.45 degrees and -23.45 degrees,
hitting zero on the equinoxes and peaking on the solstices.
'''
return 23.45 * math.sin((2 * math.pi / 365.0) * (day - 81))
In [10]:
def GetHourAngle(utc_datetime, longitude_deg):
solar_time = GetSolarTime(longitude_deg, utc_datetime)
return 15 * (12 - solar_time)
In [11]:
def GetSolarTime(longitude_deg, utc_datetime):
day = GetDayOfYear(utc_datetime)
return (((utc_datetime.hour * 60) + utc_datetime.minute + (4 * longitude_deg) + EquationOfTime(day))/60)
In [12]:
def EquationOfTime(day):
b = (2 * math.pi / 364.0) * (day - 81)
return (9.87 * math.sin(2 *b)) - (7.53 * math.cos(b)) - (1.5 * math.sin(b))
In [13]:
def get_synthetic_fsd_pri(ds):
"""
Purpose:
Calculates a time series of synthetic downwelling shortwave radiation. The
solar altitude is also output.
Useage:
qcts.get_synthetic_fsd(ds)
Author: PRI
Date: Sometime in 2014
"""
# get the latitude and longitude
lat = float(ds.globalattributes["latitude"])
lon = float(ds.globalattributes["longitude"])
# get the UTC time from the local time
ldt_UTC = qcutils.get_UTCfromlocaltime(ds)
# get the solar altitude
alt_solar = [GetAltitudeFast(lat,lon,dt) for dt in ldt_UTC]
# get the synthetic downwelling shortwave radiation
Fsd_syn = [GetRadiationDirect(dt,alt) for dt,alt in zip(ldt_UTC,alt_solar)]
Fsd_syn = numpy.ma.array(Fsd_syn)
# get the QC flag
nRecs = len(Fsd_syn)
flag = numpy.zeros(nRecs,dtype=numpy.int32)
# add the synthetic downwelling shortwave radiation to the data structure
attr = qcutils.MakeAttributeDictionary(long_name='Synthetic downwelling shortwave radiation',units='W/m2',
standard_name='surface_downwelling_shortwave_flux_in_air')
qcutils.CreateSeries(ds,"Fsd_syn",Fsd_syn,Flag=flag,Attr=attr)
# add the solar altitude to the data structure
attr = qcutils.MakeAttributeDictionary(long_name='Solar altitude',units='deg',
standard_name='not defined')
qcutils.CreateSeries(ds,"solar_altitude",alt_solar,Flag=flag,Attr=attr)
In [14]:
nc_name = qcio.get_filename_dialog()
In [15]:
ds = qcio.nc_read_series(nc_name)
In [26]:
Fsd_syn_IM = get_synthetic_fsd_im(ds)
print len(Fsd_syn_IM)
In [29]:
#Fsd_syn_PRI = get_synthetic_fsd_pri(ds)
lat = float(ds.globalattributes["latitude"])
lon = float(ds.globalattributes["longitude"])
ldt_UTC = qcutils.get_UTCfromlocaltime(ds)
alt_solar = [GetAltitudeFast(lat,lon,dt) for dt in ldt_UTC]
Fsd_syn_pri = [GetRadiationDirect(dt,alt) for dt,alt in zip(ldt_UTC,alt_solar)]
In [18]:
ldt = ds.series["DateTime"]["Data"]
Fsd,f,a = qcutils.GetSeriesasMA(ds,"Fsd")
Fsd_syn,f,a = qcutils.GetSeriesasMA(ds,"Fsd_syn")
In [30]:
fig = plt.figure()
plt.plot(ldt,Fsd,'b.')
plt.plot(ldt,Fsd_syn_pri,'r+')
plt.show()
In [27]:
print len(ldt),len(Fsd_syn_IM.flatten(order="C"))
In [20]:
fig = plt.figure()
plt.plot(ldt,Fsd,'b.')
plt.plot(ldt,Fsd_syn_IM.flatten(order="C"),'g.')
plt.show()
In [36]:
print GetOpticalDepth(300)
In [ ]: