#!/opt/EPD/bin/python
import numpy as np
import matplotlib.pyplot as plt
from optparse import OptionParser
# we are using the "option parser" to bring in optional command line arguments into our script
# for example, we can now easily change the location of the analysis to Frankfort by doing
# ./rwis.py -p fft
#
# or change the number of minutes to the nearest observation cutoff for determining valid observations to 60 min by doing
# ./rwis.py -c 60
#
# or both
# ./rwis.py -p gas -c 30
# define options and set up default values
parser = OptionParser()
parser.add_option("-p", "--prefix", dest="prefix", default="aty",
help="start of csv filenames", metavar="file name prefix")
parser.add_option("-c", "--cutoff", dest="cutoff", type="int",default=120,
help="limit of minutes for interp", metavar="minutes")
parser.add_option("-r", "--range", dest="mrange", type="int", default=30,
help="number of days", metavar="DAYS")
parser.add_option("-y", "--year", dest="myear", type="int", default=2012,
help="year", metavar="YEAR")
parser.add_option("-m", "--month", dest="mmonth", type="int", default=11,
help="month", metavar="MONTH")
# obtain information from the command line options
(options, args) = parser.parse_args()
year=options.myear
month=options.mmonth
maxday=options.mrange
pref=options.prefix
numv=options.cutoff
# set up csv file names, read data
afile = pref+'ahr.csv'
bfile = pref+'bhr.csv'
wrffile = pref+'wrf.csv'
print 'reading ',afile
print 'reading ',bfile
print 'reading ',wrffile
adata=np.genfromtxt(afile, delimiter=',')
bdata=np.genfromtxt(bfile, delimiter=',')
wrf=np.genfromtxt(wrffile, delimiter=',')
print ' total number of valid hours A site'
print adata[adata[:,5]<=numv,5].size
print ' total number of valid hours B site'
print bdata[bdata[:,5]<=numv,5].size
# compute summary error statistics
fcst=wrf[adata[:,5]<=numv,1]
obs=adata[adata[:,5]<=numv,1]
error=fcst-obs
mse=np.mean(error**2)
bias=np.mean(error)
mae=np.mean(np.abs(error))
print ' sftemp A site '
print ' mse = ',mse
print ' mae = ',mae
print ' bias = ', bias
fcst=wrf[bdata[:,5]<=numv,1]
obs=bdata[bdata[:,5]<=numv,1]
error=fcst-obs
mse=np.mean(error**2)
bias=np.mean(error)
mae=np.mean(np.abs(error))
print ' sftemp B site '
print ' mse = ',mse
print ' mae = ',mae
print ' bias = ', bias
fcst=wrf[adata[:,5]<=numv,2]
obs=adata[adata[:,5]<=numv,2]
error=fcst-obs
mse=np.mean(error**2)
bias=np.mean(error)
mae=np.mean(np.abs(error))
print ' airtemp '
print ' mse = ',mse
print ' mae = ',mae
print ' bias = ', bias
fcst=wrf[adata[:,5]<=numv,3]
obs=adata[adata[:,5]<=numv,3]
error=fcst-obs
mse=np.mean(error**2)
bias=np.mean(error)
mae=np.mean(np.abs(error))
print ' dewpoint '
print ' mse = ',mse
print ' mae = ',mae
print ' bias = ', bias
# determine the hourly errors
length = len(adata)
hr=np.zeros(24)
msehr=np.zeros(24)
maehr=np.zeros(24)
biashr=np.zeros(24)
sizehr=np.zeros(24)
mseair=np.zeros(24)
maeair=np.zeros(24)
biasair=np.zeros(24)
msedew=np.zeros(24)
maedew=np.zeros(24)
biasdew=np.zeros(24)
mseb=np.zeros(24)
maeb=np.zeros(24)
biasb=np.zeros(24)
for i in np.arange(24):
hr[i]=i+7
asubset=adata[i:length:24,:]
bsubset=bdata[i:length:24,:]
wrfsubset=wrf[i:length:24,:]
fcst=wrfsubset[asubset[:,5]<=numv,1]
obs=asubset[asubset[:,5]<=numv,1]
bobs=bsubset[asubset[:,5]<=numv,1]
airfcst=wrfsubset[asubset[:,5]<=numv,2]
airobs=asubset[asubset[:,5]<=numv,2]
dewfcst=wrfsubset[asubset[:,5]<=numv,3]
dewobs=asubset[asubset[:,5]<=numv,3]
error=fcst-obs
berror=fcst-bobs
airerror=airfcst-airobs
dewerror=dewfcst-dewobs
sizehr[i]=obs.size
msehr[i]=np.mean(error**2)
biashr[i]=np.mean(error)
maehr[i]=np.mean(np.abs(error))
mseair[i]=np.mean(airerror**2)
biasair[i]=np.mean(airerror)
maeair[i]=np.mean(np.abs(airerror))
msedew[i]=np.mean(dewerror**2)
biasdew[i]=np.mean(dewerror)
maedew[i]=np.mean(np.abs(dewerror))
mseb[i]=np.mean(berror**2)
biasb[i]=np.mean(berror)
maeb[i]=np.mean(np.abs(berror))
# plot the hourly errors
plt.figure()
plt.plot(hr,np.sqrt(msehr),'r',hr,biashr,'b',hr,maehr,'g')
plt.hold()
plt.xlabel('forecast hour')
plt.ylabel('deg C')
plt.title('A site surface temp Bias (blue) RMSE (red) MAE (green)')
plt.show()
plt.figure()
plt.plot(hr,np.sqrt(mseair),'r',hr,biasair,'b',hr,maeair,'g')
plt.hold()
plt.xlabel('forecast hour')
plt.ylabel('deg C')
plt.title('air temp Bias (blue) RMSE (red) MAE (green)')
plt.show()
plt.figure()
plt.plot(hr,np.sqrt(msedew),'r',hr,biasdew,'b',hr,maedew,'g')
plt.hold()
plt.xlabel('forecast hour')
plt.ylabel('deg C')
plt.title('dewpoint Bias (blue) RMSE (red) MAE (green)')
plt.show()
plt.figure()
plt.plot(hr,np.sqrt(mseb),'r',hr,biasb,'b',hr,maeb,'g')
plt.hold()
plt.xlabel('forecast hour')
plt.ylabel('deg C')
plt.title('B site surface temp Bias (blue) RMSE (red) MAE (green)')
plt.show()
# apply the bias found previously to the hourly forecasts, see how that improves those forecasts
mseadj=np.zeros(24)
maeadj=np.zeros(24)
biasadj=np.zeros(24)
mseairadj=np.zeros(24)
maeairadj=np.zeros(24)
biasairadj=np.zeros(24)
msedewadj=np.zeros(24)
maedewadj=np.zeros(24)
biasdewadj=np.zeros(24)
msebadj=np.zeros(24)
maebadj=np.zeros(24)
biasbadj=np.zeros(24)
for i in np.arange(24):
asubset=adata[i:length:24,:]
bsubset=bdata[i:length:24,:]
wrfsubset=wrf[i:length:24,:]
fcst=wrfsubset[asubset[:,5]<=numv,1]-biashr[i]
obs=asubset[asubset[:,5]<=numv,1]
bobs=bsubset[asubset[:,5]<=numv,1]
airfcst=wrfsubset[asubset[:,5]<=numv,2]-biasair[i]
airobs=asubset[asubset[:,5]<=numv,2]
dewfcst=wrfsubset[asubset[:,5]<=numv,3]-biasdew[i]
dewobs=asubset[asubset[:,5]<=numv,3]
error=fcst-obs
berror=fcst-bobs
airerror=airfcst-airobs
dewerror=dewfcst-dewobs
mseadj[i]=np.mean(error**2)
biasadj[i]=np.mean(error)
maeadj[i]=np.mean(np.abs(error))
mseairadj[i]=np.mean(airerror**2)
biasairadj[i]=np.mean(airerror)
maeairadj[i]=np.mean(np.abs(airerror))
msedewadj[i]=np.mean(dewerror**2)
biasdewadj[i]=np.mean(dewerror)
maedewadj[i]=np.mean(np.abs(dewerror))
msebadj[i]=np.mean(berror**2)
biasbadj[i]=np.mean(berror)
maebadj[i]=np.mean(np.abs(berror))
# plot MAE of un-adjusted (raw) forecast against MAE of bias-adjusted forecast
plt.figure()
plt.plot(hr,maehr,'b',hr,maeadj,'g')
plt.hold()
plt.xlabel('forecast hour')
plt.ylabel('deg C')
plt.title('A site raw MAE (blue) MAE of bias-adjusted fcst (green)')
plt.show()
plt.figure()
plt.plot(hr,maeair,'b',hr,maeairadj,'g')
plt.hold()
plt.xlabel('forecast hour')
plt.ylabel('deg C')
plt.title('air temp raw MAE (blue) MAE of bias-adjusted fcst (green)')
plt.show()
plt.figure()
plt.plot(hr,maedew,'b',hr,maedewadj,'g')
plt.hold()
plt.xlabel('forecast hour')
plt.ylabel('deg C')
plt.title('dewpoint raw MAE (blue) MAE of bias-adjusted fcst (green)')
plt.show()
plt.figure()
plt.plot(hr,maeb,'b',hr,maebadj,'g')
plt.hold()
plt.xlabel('forecast hour')
plt.ylabel('deg C')
plt.title('B site raw MAE (blue) MAE of bias-adjusted fcst (green)')
plt.show()