In [1]:
%run basics
%matplotlib
import pysolar
import pytz
In [2]:
site_name = "HowardSprings"
tower_name = "../Sites/"+site_name+"/Data/Portal/"+site_name+"_2014_L3.nc"
ds_tower = qcio.nc_read_series(tower_name)
site_timezone = ds_tower.globalattributes["time_zone"]
site_latitude = float(ds_tower.globalattributes["latitude"])
site_longitude = float(ds_tower.globalattributes["longitude"])
tower_timestep = int(ds_tower.globalattributes["time_step"])
print site_latitude,site_longitude
dt_tower = ds_tower.series["DateTime"]["Data"]
Fsd_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Fsd")
In [3]:
erai_name = "../ECMWF/ECMWF_2014.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)
site_lat_index = int(((latitude[0]-site_latitude)/lat_resolution)+0.5)
site_lon_index = int(((site_longitude-longitude[0])/lon_resolution)+0.5)
erai_timestep = 180
print latitude[site_lat_index],longitude[site_lon_index]
In [4]:
# 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])
#print dt_ecmwf[0],dt_ecmwf[-1]
In [5]:
# variables are dimensioned [time,latitude,longitude]
Fsd_3d = erai_file.variables["ssrd"][:,:,:]
Fsd_accum = Fsd_3d[:,site_lat_index,site_lon_index]
In [6]:
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
erai_offset = datetime.timedelta(minutes=float(erai_timestep)/2)
dt_erai_loc = [x.astimezone(site_tz) for x in dt_erai_utc]
#dt_erai_loc = [x.astimezone(site_tz) for x in dt_erai_utc]
# NOTE: will have to disable daylight saving at some stage, towers stay on Standard Time
# PRI hopes that the following line will do this ...
#dt_ecmwf_loc = [x-x.dst() for x in dt_ecmwf_loc]
# make local time timezone naive to match datetimes in OzFluxQC
dt_erai_loc_ntz = [x.replace(tzinfo=None) for x in dt_erai_loc]
dt_erai_loc_cor = [x - erai_offset for x in dt_erai_loc_ntz]
dt_erai_utc_ntz = [x.replace(tzinfo=None) for x in dt_erai_utc]
dt_erai_utc_cor = [x - erai_offset for x in dt_erai_utc_ntz]
In [7]:
# Downwelling shortwave in ERA-I is a cummulative value that is reset to 0
# at 0300 and 1500 UTC. Here we convert the cummulative values to
# 3 hourly values.
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 [8]:
fig=plt.figure()
plt.plot(dt_tower,Fsd_tower,'b-')
plt.plot(dt_erai_loc_cor,Fsd_erai_3hr,'r+')
#plt.plot(dt_erai_loc,1000*numpy.sin(numpy.deg2rad(alt_solar_3hr)),'g^')
plt.show()
In [9]:
# Now we interpolate from the 3 hourly ERA-I time step to the tower time step using the solar altitude.
nth = erai_timestep/tower_timestep
Fsd_erai_1hr = numpy.zeros(len(Fsd_erai_3hr)*3)
idx = numpy.array(range(0,len(Fsd_erai_3hr)))
for i in range(0,3):
Fsd_erai_1hr[idx*3+i] = Fsd_erai_3hr[idx]
#Fsd_erai_1hr[idx*3+1] = Fsd_erai_3hr[idx]
#Fsd_erai_1hr[idx*3+2] = Fsd_erai_3hr[idx]
# now get the 1 hour datetime series
tdhr = datetime.timedelta(minutes=60)
start_date = dt_erai_loc_cor[0]-tdhr
end_date = dt_erai_loc_cor[-1]+tdhr
dt_erai_loc_1hr = [x for x in qcutils.perdelta(start_date,end_date,tdhr)]
start_date = dt_erai_utc_cor[0]-tdhr
end_date = dt_erai_utc_cor[-1]+tdhr
dt_erai_utc_1hr = [x for x in qcutils.perdelta(start_date,end_date,tdhr)]
#dt_erai_loc_1hr = [dt_erai_loc_cor[0]-datetime.timedelta(minutes=60)]+dt_erai_loc_1hr
#dt_erai_loc_1hr = dt_erai_loc_1hr+[dt_erai_loc_cor[-1]+datetime.timedelta(minutes=60)]
In [10]:
fig=plt.figure()
plt.plot(dt_tower,Fsd_tower,'b-')
plt.plot(dt_erai_loc_1hr,Fsd_erai_1hr,'r+')
plt.show()
In [11]:
# get the solar altitude, we will use this later to interpolate the ERA Interim
# data from the ERA-I 3 hour time step to the tower time step.
# alt_solar is in degrees
alt_solar_3hr = numpy.array([pysolar.GetAltitude(site_latitude,site_longitude,dt) for dt in dt_erai_utc_cor])
In [12]:
#alt_solar_3hr = numpy.ma.masked_where(alt_solar_3hr<0,alt_solar_3hr)
coef_3hr = Fsd_erai_3hr/numpy.sin(numpy.deg2rad(alt_solar_3hr))
In [13]:
fig = plt.figure()
ax1 = plt.subplot(311)
plt.plot(dt_erai_loc_cor,alt_solar_3hr,'bo')
ax2 = plt.subplot(312,sharex=ax1)
plt.plot(dt_erai_loc_cor,Fsd_erai_3hr,'bo')
ax3 = plt.subplot(313,sharex=ax1)
plt.plot(dt_erai_loc_cor,coef_3hr,'bo')
plt.show()
In [17]:
from scipy.interpolate import InterpolatedUnivariateSpline
alt_solar_1hr = numpy.array([pysolar.GetAltitude(site_latitude,site_longitude,dt) for dt in dt_erai_utc_1hr])
idx = numpy.where(alt_solar_1hr<=0)[0]
alt_solar_1hr[idx] = float(0)
#coef_1hr = numpy.ma.zeros(len(coef_3hr)*3)
#idx = numpy.array(range(0,len(coef_3hr)))
#coef_1hr[idx*3] = coef_3hr[idx]
#coef_1hr[idx*3+1] = coef_3hr[idx]
#coef_1hr[idx*3+2] = coef_3hr[idx]
idx_3hr = numpy.arange(0,len(coef_3hr))*3
idx_1hr = numpy.arange(0,len(coef_3hr)*3)
s = InterpolatedUnivariateSpline(idx_3hr, coef_3hr, k=1)
coef_1hr = s(idx_1hr)
Fsd_erai_1hr = coef_1hr*numpy.sin(numpy.deg2rad(alt_solar_1hr))
In [18]:
fig = plt.figure()
plt.plot(dt_tower,Fsd_tower,'b-')
plt.plot(dt_erai_loc_cor,Fsd_erai_3hr,'r+')
plt.plot(dt_erai_loc_1hr,Fsd_erai_1hr,'g^')
plt.tight_layout()
plt.show()
In [ ]:
# now get the true solar altitude
start_date = dt_erai_utc_cor[0]
end_date = dt_erai_utc_cor[-1]
dt_erai_utc_1hr = [x for x in qcutils.perdelta(start_date,end_date,tdhr)]
dt_erai_utc_1hr = [dt_erai[0]-datetime.timedelta(minutes=60)]+dt_erai_utc_1hr
dt_erai_utc_1hr = dt_erai_utc_1hr+[dt_erai[-1]+datetime.timedelta(minutes=60)]
alt_solar = numpy.array([pysolar.GetAltitude(site_latitude,site_longitude,dt) for dt in dt_erai_utc_1hr])
In [ ]:
fig=plt.figure()
plt.plot(dt_erai_cor,alt_solar_3hr)
plt.plot(dt_erai_1hr,alt_solar_1hr,'r+')
plt.plot(dt_erai_1hr,alt_solar,'g^')
plt.show()
In [ ]:
fig=plt.figure()
plt.plot(dt_erai_cor,1000*numpy.sin(numpy.deg2rad(alt_solar_3hr)),'b.')
#plt.plot(dt_erai_1hr,1000*numpy.sin(numpy.deg2rad(alt_solar)),'b^')
plt.plot(dt_erai_cor,Fsd_erai_3hr,'r+')
plt.show()
In [ ]:
ratio = numpy.deg2rad(alt_solar)/numpy.deg2rad(alt_solar_1hr)
Fsd_erai_1hr = ratio*Fsd_erai_1hr
fig=plt.figure()
plt.plot(dt_tower,Fsd_tower,'b-')
plt.plot(dt_erai_1hr,Fsd_erai_1hr,'r+')
plt.ylim([0,1200])
plt.show()
In [ ]:
from scipy.interpolate import InterpolatedUnivariateSpline
order = 1
s = InterpolatedUnivariateSpline(xi, yi, k=order)
y = s(x)