In [1]:
%run basics
%matplotlib
import pysolar
import pytz
from scipy.interpolate import InterpolatedUnivariateSpline


Using matplotlib backend: Qt4Agg

In [2]:
do_plot = True

In [61]:
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"])
dt_tower = ds_tower.series["DateTime"]["Data"]
Fsd_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Fsd")
Fsu_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Fsu")
Fn_sw_tower = Fsd_tower - Fsu_tower
Fld_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Fld")
Flu_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Flu")
Fn_lw_tower = Fld_tower - Flu_tower
Fn_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Fn")
Fh_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Fh")
Fe_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Fe")
Fg_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Fg")
Fa_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Fa")
ps_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"ps")
Ta_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Ta")
Ah_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Ah")
Precip_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Precip")
Sws_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Sws")
Ts_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Ts")
Ws_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Ws")
Wd_tower,f,a = qcutils.GetSeriesasMA(ds_tower,"Wd")

In [4]:
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)
# index of the site in latitude dimension
site_lat_index = int(((latitude[0]-site_latitude)/lat_resolution)+0.5)
erai_latitude = latitude[site_lat_index]
# index of the site in longitude dimension
site_lon_index = int(((site_longitude-longitude[0])/lon_resolution)+0.5)
erai_longitude = longitude[site_lon_index]
erai_timestep = 180

In [5]:
# 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])
# 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.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_utc_ntz = [x.replace(tzinfo=None) for x in dt_erai_utc]
# get the datetime in the middle of the accumulation period
erai_offset = datetime.timedelta(minutes=float(erai_timestep)/2)
# subtract this from the ERA-I datetime to align ERA-I 3 hourly points with the middle
# of the accumulation period
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]
# now we get the datetime series at the tower time step
tdts = datetime.timedelta(minutes=tower_timestep)
# local time for plotting
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)]
# UTC for calculating solar altitude
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)]
# UTC at tower time step as number for interpolation
erai_time_tts = netCDF4.date2num(dt_erai_utc_tts,time_units)
erai_time_3hr = netCDF4.date2num(dt_erai_utc_cor,time_units)
# 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.
# NOTE: alt_solar is in degrees
alt_solar_3hr = numpy.array([pysolar.GetAltitude(erai_latitude,erai_longitude,dt) for dt in dt_erai_utc_cor])
# get the solar altitude at the tower time step
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 [6]:
# Interpolate the 3 hourly accumulated downwelling shortwave to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Fsd_3d = erai_file.variables["ssrd"][:,:,:]
Fsd_accum = Fsd_3d[:,site_lat_index,site_lon_index]
# 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)
# normalise the ERA-I downwelling shortwave by the solar altitude
coef_3hr = Fsd_erai_3hr/numpy.sin(numpy.deg2rad(alt_solar_3hr))
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, coef_3hr, k=1)
# get the coefficient at the tower time step
coef_tts = s(erai_time_tts)
# get the downwelling solar radiation at the tower time step
Fsd_erai_tts = coef_tts*numpy.sin(numpy.deg2rad(alt_solar_tts))
if do_plot==True:
    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_tts,Fsd_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [15]:
# Interpolate the 3 hourly accumulated net shortwave to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Fn_sw_3d = erai_file.variables["ssr"][:,:,:]
Fn_sw_accum = Fn_sw_3d[:,site_lat_index,site_lon_index]
# Net 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.
Fn_sw_erai_3hr = numpy.ediff1d(Fn_sw_accum,to_begin=0)
# deal with the reset times at 0300 and 1500
idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]
Fn_sw_erai_3hr[idx] = Fn_sw_accum[idx]
# get the average value over the 3 hourly period
Fn_sw_erai_3hr = Fn_sw_erai_3hr/(erai_timestep*60)
# normalise the ERA-I et shortwave by the solar altitude
coef_3hr = Fn_sw_erai_3hr/numpy.sin(numpy.deg2rad(alt_solar_3hr))
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, coef_3hr, k=1)
# get the coefficient at the tower time step
coef_tts = s(erai_time_tts)
# get the downwelling solar radiation at the tower time step
Fn_sw_erai_tts = coef_tts*numpy.sin(numpy.deg2rad(alt_solar_tts))
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Fn_sw_tower,'b-')
    plt.plot(dt_erai_loc_cor,Fn_sw_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,Fn_sw_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [16]:
Fsu_erai_tts = Fsd_erai_tts - Fn_sw_erai_tts
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Fsu_tower,'b-')
    plt.plot(dt_erai_loc_tts,Fsu_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [17]:
# Interpolate the 3 hourly accumulated downwelling longwave to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Fld_3d = erai_file.variables["strd"][:,:,:]
Fld_accum = Fld_3d[:,site_lat_index,site_lon_index]
# Downwelling longwave 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.
Fld_erai_3hr = numpy.ediff1d(Fld_accum,to_begin=0)
idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]
Fld_erai_3hr[idx] = Fld_accum[idx]
Fld_erai_3hr = Fld_erai_3hr/(erai_timestep*60)
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Fld_erai_3hr, k=1)
# get the downwelling longwave at the tower time step
Fld_erai_tts = s(erai_time_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Fld_tower,'b-')
    plt.plot(dt_erai_loc_cor,Fld_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,Fld_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [18]:
# Interpolate the 3 hourly accumulated net longwave to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Fn_lw_3d = erai_file.variables["str"][:,:,:]
Fn_lw_accum = Fn_lw_3d[:,site_lat_index,site_lon_index]
# Net longwave in ERA-I is a cummulative value that is reset at 0300 and 1500 UTC.
# Here we convert the cummulative values to 3 hourly values.
Fn_lw_erai_3hr = numpy.ediff1d(Fn_lw_accum,to_begin=0)
# deal with the reset times at 0300 and 1500
idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]
Fn_lw_erai_3hr[idx] = Fn_lw_accum[idx]
# get the average value over the 3 hourly period
Fn_lw_erai_3hr = Fn_lw_erai_3hr/(erai_timestep*60)
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Fn_lw_erai_3hr, k=1)
# get the net longwave at the tower time step
Fn_lw_erai_tts = s(erai_time_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Fn_lw_tower,'b-')
    plt.plot(dt_erai_loc_cor,Fn_lw_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,Fn_lw_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [19]:
Flu_erai_tts = Fld_erai_tts - Fn_lw_erai_tts
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Flu_tower,'b-')
    plt.plot(dt_erai_loc_tts,Flu_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [20]:
Fn_erai_tts = (Fsd_erai_tts - Fsu_erai_tts) + (Fld_erai_tts - Flu_erai_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Fn_tower,'b-')
    plt.plot(dt_erai_loc_tts,Fn_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [64]:
# Interpolate the 3 hourly accumulated sensible heat flux to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Fh_3d = erai_file.variables["sshf"][:,:,:]
Fh_accum = float(-1)*Fh_3d[:,site_lat_index,site_lon_index]
# Sensible heat flux in ERA-I is a cummulative value that is reset at 0300 and 1500 UTC.
# Here we convert the cummulative values to 3 hourly values.
Fh_erai_3hr = numpy.ediff1d(Fh_accum,to_begin=0)
# deal with the reset times at 0300 and 1500
idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]
Fh_erai_3hr[idx] = Fh_accum[idx]
# get the average value over the 3 hourly period
Fh_erai_3hr = Fh_erai_3hr/(erai_timestep*60)
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Fh_erai_3hr, k=1)
# get the net longwave at the tower time step
Fh_erai_tts = s(erai_time_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Fh_tower,'b-')
    plt.plot(dt_erai_loc_cor,Fh_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,Fh_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [65]:
# Interpolate the 3 hourly accumulated latent heat flux to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Fe_3d = erai_file.variables["slhf"][:,:,:]
Fe_accum = float(-1)*Fe_3d[:,site_lat_index,site_lon_index]
# Latent heat flux in ERA-I is a cummulative value that is reset at 0300 and 1500 UTC.
# Here we convert the cummulative values to 3 hourly values.
Fe_erai_3hr = numpy.ediff1d(Fe_accum,to_begin=0)
# deal with the reset times at 0300 and 1500
idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]
Fe_erai_3hr[idx] = Fe_accum[idx]
# get the average value over the 3 hourly period
Fe_erai_3hr = Fe_erai_3hr/(erai_timestep*60)
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Fe_erai_3hr, k=1)
# get the net longwave at the tower time step
Fe_erai_tts = s(erai_time_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Fe_tower,'b-')
    plt.plot(dt_erai_loc_cor,Fe_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,Fe_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [68]:
# get Fg as residual
Fg_erai_tts = Fn_erai_tts - Fh_erai_tts - Fe_erai_tts
# and then Fa
Fa_erai_tts = Fn_erai_tts - Fg_erai_tts
if do_plot==True:
    fig = plt.figure()
    ax1 = plt.subplot(211)
    plt.plot(dt_tower,Fg_tower,'b-')
    plt.plot(dt_erai_loc_tts,Fg_erai_tts,'g^')
    ax2 = plt.subplot(212,sharex=ax1)
    plt.plot(dt_tower,Fa_tower,'b-')
    plt.plot(dt_erai_loc_tts,Fa_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [25]:
# Interpolate the 3 hourly air pressure to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
ps_3d = erai_file.variables["sp"][:,:,:]
ps_erai_3hr = ps_3d[:,site_lat_index,site_lon_index]/float(1000)
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, ps_erai_3hr, k=1)
# get the air pressure at the tower time step
ps_erai_tts = s(erai_time_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,ps_tower,'b-')
    plt.plot(dt_erai_loc_cor,ps_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,ps_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [22]:
# Interpolate the 3 hourly air temperature to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Ta_3d = erai_file.variables["t2m"][:,:,:]
Ta_erai_3hr = Ta_3d[:,site_lat_index,site_lon_index] - 273.15
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Ta_erai_3hr, k=1)
# get the air temperature at the tower time step
Ta_erai_tts = s(erai_time_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Ta_tower,'b-')
    plt.plot(dt_erai_loc_cor,Ta_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,Ta_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [26]:
# Interpolate the 3 hourly dew point temperature to the tower time step
# and convert to Ah, RH and q
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Td_3d = erai_file.variables["d2m"][:,:,:]
Td_erai_3hr = Td_3d[:,site_lat_index,site_lon_index] - 273.15
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Td_erai_3hr, k=1)
# get the dew point temperature at the tower time step
Td_erai_tts = s(erai_time_tts)
# get the relative humidity
es_erai_tts = mf.es(Ta_erai_tts)
e_erai_tts = mf.es(Td_erai_tts)
RH_erai_tts = float(100)*e_erai_tts/es_erai_tts
# get the absolute humidity
Ah_erai_tts = mf.absolutehumidityfromRH(Ta_erai_tts,RH_erai_tts)
# get the specific humidity
q_erai_tts = mf.specifichumidityfromRH(RH_erai_tts,Ta_erai_tts,ps_erai_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Ah_tower,'b-')
    plt.plot(dt_erai_loc_tts,Ah_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [27]:
# Interpolate the 3 hourly boundary layer height to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Habl_3d = erai_file.variables["blh"][:,:,:]
Habl_erai_3hr = Habl_3d[:,site_lat_index,site_lon_index]
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Habl_erai_3hr, k=1)
# get the boundary layer height at the tower time step
Habl_erai_tts = s(erai_time_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_erai_loc_cor,Habl_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,Habl_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [51]:
# Spread the 3 hourly accumulated precipitation to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Precip_3d = erai_file.variables["tp"][:,:,:]
Precip_accum = Precip_3d[:,site_lat_index,site_lon_index]

Precip_erai_3hr = numpy.ediff1d(Precip_accum,to_begin=0)
idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]
Precip_erai_3hr[idx] = Precip_accum[idx]
Precip_erai_3hr = Precip_erai_3hr*float(1000)

Precip_erai_tts = numpy.zeros(len(dt_erai_loc_tts))
idx = qcutils.FindIndicesOfBInA(dt_erai_loc_cor,dt_erai_loc_tts)
Precip_erai_tts[idx] = Precip_erai_3hr

if do_plot:
    fig = plt.figure()
    plt.plot(dt_tower,Precip_tower,'b-')
    plt.plot(dt_erai_loc_cor,Precip_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,Precip_erai_tts,'g^')
    plt.show()

#idx = qcutils.FindIndicesOfBInA(dt_erai_loc_tts,dt_tower)

#Precip_erai_tts_ma = numpy.ma.array(Precip_erai_tts,mask=Precip_tower.mask)
#Precip_erai_tts_cs = numpy.ma.cumsum(Precip_erai_tts)
#Precip_tower_cs = numpy.ma.cumsum(Precip_tower)

#if do_plot:
#    fig = plt.figure()
#    plt.plot(dt_tower,Precip_tower_cs,'b-')
#    plt.plot(dt_erai_loc_tts,Precip_erai_tts_cs,'g^')
#    plt.show()

In [53]:
# Interpolate the 3 hourly soil moisture to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Sws_3d = erai_file.variables["swvl1"][:,:,:]
Sws_erai_3hr = Sws_3d[:,site_lat_index,site_lon_index]
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Sws_erai_3hr, k=1)
# get the soil moisture at the tower time step
Sws_erai_tts = s(erai_time_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Sws_tower,'b-')
    plt.plot(dt_erai_loc_cor,Sws_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,Sws_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [57]:
# Interpolate the 3 hourly soil temperature to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Ts_3d = erai_file.variables["stl1"][:,:,:]
Ts_erai_3hr = Ts_3d[:,site_lat_index,site_lon_index] - 273.15
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Ts_erai_3hr, k=1)
# get the soil moisture at the tower time step
Ts_erai_tts = s(erai_time_tts)
if do_plot==True:
    fig = plt.figure()
    plt.plot(dt_tower,Ts_tower,'b-')
    plt.plot(dt_erai_loc_cor,Ts_erai_3hr,'r+')
    plt.plot(dt_erai_loc_tts,Ts_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [60]:
# Interpolate the 3 hourly U and V components to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
# U first ...
U_3d = erai_file.variables["u10"][:,:,:]
U_erai_3hr = U_3d[:,site_lat_index,site_lon_index]
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, U_erai_3hr, k=1)
# get the soil moisture at the tower time step
U_erai_tts = s(erai_time_tts)
# ... then V
V_3d = erai_file.variables["v10"][:,:,:]
V_erai_3hr = V_3d[:,site_lat_index,site_lon_index]
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, V_erai_3hr, k=1)
# get the soil moisture at the tower time step
V_erai_tts = s(erai_time_tts)
# now get the wind speed and direction
Ws_erai_tts = numpy.sqrt(U_erai_tts*U_erai_tts + V_erai_tts*V_erai_tts)
Wd_erai_tts = float(270) - numpy.arctan2(V_erai_tts,U_erai_tts)*float(180)/numpy.pi
idx = numpy.where(Wd_erai_tts>360)[0]
if len(idx)>0: Wd_erai_tts[idx] = Wd_erai_tts[idx] - float(360)

if do_plot==True:
    fig = plt.figure()
    ax1 = plt.subplot(211)
    plt.plot(dt_tower,Ws_tower,'b-')
    plt.plot(dt_erai_loc_tts,Ws_erai_tts,'g^')
    ax2 = plt.subplot(212,sharex=ax1)
    plt.plot(dt_tower,Wd_tower,'b-')
    plt.plot(dt_erai_loc_tts,Wd_erai_tts,'g^')
    plt.tight_layout()
    plt.show()

In [ ]: