Skip to content

GitLab

Commit 575b127c authored by saeed's avatar saeed
Browse files

improvement on sshval and taylordiagram

parent 08f764e3
Hide whitespace changes
Inline Side-by-side
sshval.py
View file @ 575b127c
......@@ -15,7 +15,7 @@ cstyle = []
ccolor = []
rms = []
cor = []
rmsd = []
def format_time_string(time_str):
"""
Converts 2017-10-01 20:00 -> 201710012000
......@@ -80,6 +80,7 @@ def readConfig():
with open(config_file, 'r') as fptr_cfg:
cfg = json.load(fptr_cfg)
server = cfg["server"]
dovalidation = cfg["dovalidation"]
station = cfg["station"]
startdate = format_time_string(cfg["startdate"])
......@@ -88,13 +89,15 @@ def readConfig():
expname = cfg["experiment"]
oper = cfg["operational"]
tickint = cfg["tickinterval_hour"]
return dovalidation, station, startdate, enddate,obstyle, expname, oper, tickint
return server,dovalidation, station, startdate, enddate,obstyle, expname, oper, tickint
def readObs(station, startdate,enddate,obstyle):
param = {'from': startdate, 'too': enddate}
res = req.get("http://oceandata.smhi.se/ssh/"+station+"/OBSERVATION",params = param)
dobs = res.json()
def readObs(server,station, startdate,enddate,obstyle):
if server["prod"]: param = {'from': startdate, 'too': enddate}
if server["utv"]:param = {'highfreq':'false','from': startdate, 'too': enddate}
if server["prod"]: res = req.get("http://oceandata.smhi.se/ssh/"+station+"/OBSERVATION",params = param)
if server["utv"]:res = req.get("http://oceandata-utv.smhi.se/ssh/"+station+"/OBSERVATION",params = param)
dobs = res.json()
date_obs = dobs.keys()
dmiss = fillObs(date_obs,startdate, enddate)
if dmiss.keys():
......@@ -107,6 +110,7 @@ def readObs(station, startdate,enddate,obstyle):
obstd = np.std(vobsm)
rms.append(round(obstd,3))
cor.append(1.0)
rmsd.append(0.0)
clegend.append("OBSERVATION")
vssh.append(vobsm)
cstyle.append(obstyle["line"])
......@@ -126,18 +130,22 @@ def readExpr(expname,station,startdate,enddate,vobsm):
lrms = np.sqrt(np.mean((vmod - vobsm)**2),dtype=np.float64)
lrms = np.sqrt(((vmod-np.mean(vmod)-vobsm+ np.mean(vobsm))**2).mean())
obstd = np.std(vobsm)
modstd = np.std(vmod - vobsm)
rms.append(round(modstd/obstd,3))
difstd = np.std(vmod - vobsm)
modstd = np.std(vmod)
rms.append(round(modstd,3))
cor.append(round(lcorr,3))
rmsd.append(round(difstd,3))
clegend.append(model)
vssh.append(vmod)
cstyle.append(expname[model]["line"])
ccolor.append(expname[model]["color"])
def readOper(oper,station,startdate,enddate,vobsm):
param = {'from': startdate, 'too': enddate}
def readOper(server,oper,station,startdate,enddate,vobsm):
if server["prod"]:param = {'from': startdate, 'too': enddate}
if server["utv"]: param = {'highfreq':'false','from': startdate, 'too': enddate}
for model in oper:
res1 = req.get("http://oceandata.smhi.se/ssh/"+station+"/"+model,params = param)
if server["prod"]:res1 = req.get("http://oceandata.smhi.se/ssh/"+station+"/"+model,params = param)
if server["utv"]:res1 = req.get("http://oceandata-utv.smhi.se/ssh/"+station+"/"+model,params = param)
doper = res1.json()
vm = pd.DataFrame.from_dict(doper,orient="index")
voper = vm["raw"].loc[startdate:enddate].values
......@@ -146,23 +154,24 @@ def readOper(oper,station,startdate,enddate,vobsm):
lrms = np.sqrt(np.mean((voper - vobsm)**2),dtype=np.float64)
lrms = np.sqrt(((voper-np.mean(voper)-vobsm+ np.mean(vobsm))**2).mean())
lrms = np.sqrt(((voper-vobsm)**2).mean())
modstd = np.std(voper - vobsm)
difstd = np.std(voper - vobsm)
modstd = np.std(voper)
obstd = np.std(vobsm)
rms.append(round(modstd/obstd,3))
rms.append(round(modstd,3))
cor.append(round(lcorr,3))
rmsd.append(round(difstd,3))
clegend.append(oper[model]["title"])
vssh.append(voper)
cstyle.append(oper[model]["line"])
ccolor.append(oper[model]["color"])
##############################
def main():
dovalidation, station, startdate, enddate, obstyle, expname, oper,tickint = readConfig()
server,dovalidation, station, startdate, enddate, obstyle, expname, oper,tickint = readConfig()
vobsm = readObs(station, startdate, enddate,obstyle)
vobsm = readObs(server,station, startdate, enddate,obstyle)
if dovalidation["experiment"] :readExpr(expname,station,startdate,enddate,vobsm)
if dovalidation["operational"] :readOper(oper,station,startdate,enddate,vobsm)
if dovalidation["operational"] :readOper(server,oper,station,startdate,enddate,vobsm)
stick = getDateRange(startdate,enddate)
itick = range(len(stick))
......@@ -187,7 +196,7 @@ def main():
plt.close(1)
TaylorDiagram(RMSVEC=rms, CORVEC=cor,COLORVEC=ccolor,LABELVEC=clegend, station=station, info=startdate+"-"+enddate)
TaylorDiagram(RMSVEC=rms, RMSDVEC=rmsd, CORVEC=cor,COLORVEC=ccolor,LABELVEC=clegend, station=station, info=startdate+"-"+enddate)
if __name__ == "__main__":
main()
taylorDiagram.py
View file @ 575b127c
def TaylorDiagram(RMSVEC, CORVEC,COLORVEC,LABELVEC, station, info):
def TaylorDiagram(RMSVEC, RMSDVEC, CORVEC,COLORVEC,LABELVEC, station, info):
import numpy as np
import matplotlib.pyplot as plt
import math
rms_max = float(math.ceil(max(RMSVEC)))
if rms_max <= 1.0: delta = 0.2
if rms_max >= 2.0 and rms_max <= 3.0: delta = 0.25
if rms_max >=4 and rms_max <= 5.0: delta = 0.5
if rms_max >=6 and rms_max <= 7.0: delta = 0.5
if rms_max >=8 and rms_max <= 9.0: delta = 1.0
if rms_max >=10: delta = 1.0
X=np.arange(0.0,rms_max+delta/100.0,delta/100.0)
Y=np.arange(0.0,rms_max+delta/100.0,delta/100.0)
print "name",LABELVEC
print "cor",CORVEC
print "rms ",RMSVEC
print "rmsd",RMSDVEC
rms_max = max(RMSVEC)
delta = rms_max/10.0
X=np.arange(0.0,rms_max+delta,delta/100.0)
Y=np.arange(0.0,rms_max+delta,delta/100.0)
nx=X.shape[0]
ny=Y.shape[0]
h = np.zeros((ny,nx),dtype=np.float32)
for j in range(ny):
for i in range(nx):
h[j,i] = math.sqrt(X[i] * X[i] + Y[j] * Y[j])
##########################################
rmsd_max = max(RMSDVEC)
ddelta = rmsd_max/10.0
XX=np.arange(0.0,rms_max+delta,delta/100.0)
YY=np.arange(0.0,rms_max+delta,delta/100.0)
nxx=XX.shape[0]
nyy=YY.shape[0]
hh = np.zeros((nyy,nxx),dtype=np.float32)
hh[:,:] = -1.0
for j in range(nyy):
for i in range(nxx):
if (math.sqrt((XX[i])*(XX[i])+ YY[j]*YY[j]) <= rms_max):hh[j,i] = math.sqrt((XX[i]-RMSVEC[0]) * (XX[i]-RMSVEC[0]) + YY[j] * YY[j])
#hh[j,i] = math.sqrt((XX[i]-RMSVEC[0]) * (XX[i]-RMSVEC[0]) + YY[j] * YY[j])
hh=np.ma.masked_where(hh==-1.0,hh)
###########################################
fig=plt.figure(num=1,figsize=(9.0,9.0),dpi=300,facecolor='w',edgecolor='k')
ax = fig.add_axes([0.08, 0.08, 0.8, 0.8], axisbg = '1.0')
ax.set_xlabel('RMSE/STD_OBS',fontsize='15',weight='bold',color="green")
ax.set_ylabel('RMSE/STD_OBS',fontsize='15',weight='bold',color="green")
ax.set_xlabel(' ',fontsize='15',weight='bold',color="green")
ax.set_ylabel('Standard Deviation',fontsize='15',weight='bold',color="green")
ax.grid(False)
vc=np.arange(0.0,rms_max+delta,delta)
ax.contour(X,Y,h,vc,colors='0.5',linewidths=0.2)
nl, = vc.shape
lines=[]
for i in range(nl-1):
lines.append("dashed")
lines.append("solid")
ax.contour(X,Y,h,vc,colors='0.5',linestyles=lines,linewidths=0.2)
vvc=np.arange(0.0,rmsd_max+ddelta,ddelta)
ax.contour(XX,YY,hh,vvc,colors='0.5',linestyles="solid",linewidths=0.3)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_xlim(0.0,rms_max + 0.005)
ax.set_ylim(0.0,rms_max + 0.005)
ax.set_xlim(0.0,rms_max+delta)
ax.set_ylim(0.0,rms_max+delta)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
radius = np.arange(0.0,rms_max + delta,delta)
xangle = list(np.arange(0.0,0.9,0.10)) + [0.9, 0.95, 0.99]
rdmax = np.amax(radius)
for rd in radius:
for ang in xangle:
ax.plot([0.0,rd*ang],[0.0,rd*math.sqrt(1.0 - (ang * ang))],color="0.5",ls="-",lw="0.10")
if rd == rms_max and ang > 0.0: ax.text(rd*ang, rd*math.sqrt(1.0 - (ang * ang)), str(ang),fontsize=10)
if rd == rdmax and ang > 0.0: ax.text(rd*ang, rd*math.sqrt(1.0 - (ang * ang)), str(ang),fontsize=10)
ang = 0.65
ax.text((1.040)*rms_max*ang, (1.04)*rms_max*math.sqrt(1.0 - (ang * ang)), "Correlation",color="Steelblue",fontsize=20,rotation=-45)
......@@ -52,9 +75,9 @@ def TaylorDiagram(RMSVEC, CORVEC,COLORVEC,LABELVEC, station, info):
for ii in range(lenm):
vrms, vcor, vcol, vlab = RMSVEC[ii], CORVEC[ii], COLORVEC[ii],LABELVEC[ii]
if ii == 0 :
#line, = ax.plot(vrms*vcor,vrms*math.sqrt(1.0 - (vcor * vcor)),'o',color=vcol, label=vlab,ms=8)
#line.set_clip_on(False)
ax.text(rms_max/10.0, -rms_max/8.0, "STD_OBS "+str(vrms),color="red",fontsize=12,rotation=0)
line, = ax.plot(vrms*vcor,vrms*math.sqrt(1.0 - (vcor * vcor)),'o',color=vcol, label=vlab,ms=8)
line.set_clip_on(False)
#ax.text(rms_max/10.0, -rms_max/8.0, "STD_OBS "+str(vrms),color="red",fontsize=12,rotation=0)
else:ax.plot(vrms*vcor,vrms*math.sqrt(1.0 - (vcor * vcor)),'o',color=vcol, label=vlab,ms=8)
ax.legend(numpoints=1,loc = 'best',prop=dict(size='small'),fontsize=12)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment