In [1]:
%run basics
%matplotlib
import pytz
import pysolar
from scipy.interpolate import InterpolatedUnivariateSpline
In [2]:
erai_name = "../ERAI/ERAI_2012.nc"
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 [3]:
# details for Whroo
site_name = "Whroo"
site_timezone = "Australia/Melbourne"
site_latitude = -36.67305
site_longitude = 145.02621
site_timestep = 30
site_lat_index = int(((latitude[0]-site_latitude)/lat_resolution)+0.5)
erai_latitude = latitude[site_lat_index]
site_lon_index = int(((site_longitude-longitude[0])/lon_resolution)+0.5)
erai_longitude = longitude[site_lon_index]
print site_latitude,erai_latitude
print site_longitude,erai_longitude
In [4]:
# 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-x.dst() for x in dt_erai_loc]
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]
erai_offset = datetime.timedelta(minutes=float(erai_timestep)/2)
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]
tdts = datetime.timedelta(minutes=site_timestep)
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)]
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)]
erai_time_tts = netCDF4.date2num(dt_erai_utc_tts,time_units)
erai_time_3hr = netCDF4.date2num(dt_erai_utc_cor,time_units)
alt_solar_3hr = numpy.array([pysolar.GetAltitude(erai_latitude,erai_longitude,dt) for dt in dt_erai_utc_cor])
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)
In [5]:
Fsd_3d = erai_file.variables["ssrd"][:,:,:]
Fsd_accum = Fsd_3d[:,site_lat_index,site_lon_index]
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)
In [ ]:
fig = plt.figure()
plt.plot(dt_erai_loc_cor,Fsd_erai_3hr)
plt.show()
In [18]:
#sa = numpy.array(alt_solar_3hr,copy=True)
#idx = numpy.where(sa<=float(5))[0]
#sa[idx] = float(5)
fives = float(5)*numpy.ones(len(alt_solar_3hr))
sa = numpy.where(alt_solar_3hr<=float(5),fives,alt_solar_3hr)
coef_3hr = Fsd_erai_3hr/numpy.sin(numpy.deg2rad(sa))
s = InterpolatedUnivariateSpline(erai_time_3hr,coef_3hr,k=1)
coef_tts = s(erai_time_tts)
Fsd_erai_tts = coef_tts*numpy.sin(numpy.deg2rad(alt_solar_tts))
In [19]:
fig=plt.figure()
ax1 = plt.subplot(311)
ax1.plot(dt_erai_loc_cor,coef_3hr,'ro')
ax1.plot(dt_erai_loc_tts,coef_tts,'b+')
ax2 = plt.subplot(312,sharex=ax1)
ax2.plot(dt_erai_loc_cor,numpy.sin(numpy.deg2rad(alt_solar_3hr)),'ro')
ax3 = plt.subplot(313,sharex=ax1)
ax3.plot(dt_erai_loc_cor,Fsd_erai_3hr,'ro')
ax3.plot(dt_erai_loc_tts,Fsd_erai_tts,'b+')
plt.show()
In [13]:
print numpy.sin(numpy.deg2rad(5))
In [ ]: