Skip to content
GitLab
Commit e216a845 authored by saeed's avatar saeed
Browse files

taylor and sshval are in the very good shape now

parent 10cf7525
Hide whitespace changes
Inline Side-by-side
sshval.py
View file @ e216a845
......@@ -5,13 +5,14 @@ import requests as req
import matplotlib.pyplot as plt
import json
import datetime
from scipy import stats
from taylorDiagram import TaylorDiagram
reload(sys)
sys.setdefaultencoding('utf8')
clegend = []
vssh = []
cstyle = []
ccolor = []
rms = []
cor = []
......@@ -20,12 +21,10 @@ def format_time_string(time_str):
Converts 2017-10-01 20:00 -> 201710012000
"""
input_format = '%Y-%m-%d %H:%M'
output_format = '%Y%m%d%H'
return datetime.datetime.strptime(time_str, input_format).strftime(output_format)
def strDate(day):
cyear, cmonth, cday, chour="%04d"%int(day.year),"%02d"%int(day.month),"%02d"%int(day.day),"%02d"%int(day.hour)
return "".join([cyear, cmonth, cday,chour])
......@@ -42,13 +41,30 @@ def getDateRange(fromt, tot):
return dateticks
def fillObs(date_obs,fromt, tot):
initd = datetime.datetime(int(fromt[0:4]),int(fromt[4:6]),int(fromt[6:8]),int(fromt[8:10]))
totd = datetime.datetime(int(tot[0:4]),int(tot[4:6]),int(tot[6:8]),int(tot[8:10]))
step = datetime.timedelta(hours=1)
dmiss = {}
dmcor = {}
while initd <= totd:
if strDate(initd) not in date_obs:
dmiss[strDate(initd)] = {}
dmcor[strDate(initd)] = {}
dmcor[strDate(initd)]["raw"] = -999.0
dmiss[strDate(initd)]["raw"] = np.nan
initd = initd + step
return dmcor
def biasMean(vobs, vmodel):
if vobs.shape[0] != vmodel.shape[0] :
print "vobs and vmodel should have the same size"
sys.exit()
return stats.nanmean(vobs- vmodel)
return np.mean(vobs- vmodel)
def readConfig():
......@@ -78,78 +94,71 @@ def readConfig():
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()
dobs = res.json()
date_obs = dobs.keys()
dmiss=fillObs(date_obs,startdate, enddate)
dmiss = fillObs(date_obs,startdate, enddate)
if dmiss.keys():
dobs.update(dmiss)
vo = pd.DataFrame.from_dict(dobs,orient="index")
vob = vo.loc[startdate:enddate].values
vobs = np.squeeze(vob)
voc = pd.DataFrame.from_dict(dobs,orient="index")
vobc = voc.loc[startdate:enddate].values
vobsc = np.squeeze(vobc)
vobsm = np.ma.masked_where(vobsc==-999.0,vobsc)
clegend.append("OBSERVATION")
vssh.append(vobs)
cstyle.append(obstyle)
vssh.append(vobsm)
cstyle.append(obstyle["line"])
ccolor.append(obstyle["color"])
return vobs
def fillObs(date_obs,fromt, tot):
initd = datetime.datetime(int(fromt[0:4]),int(fromt[4:6]),int(fromt[6:8]),int(fromt[8:10]))
totd = datetime.datetime(int(tot[0:4]),int(tot[4:6]),int(tot[6:8]),int(tot[8:10]))
step = datetime.timedelta(hours=1)
dmiss = {}
while initd <= totd:
if strDate(initd) not in date_obs:
dmiss[strDate(initd)] = {}
dmiss[strDate(initd)]["raw"] = np.nan
initd = initd + step
return dmiss
return vobsm
def readExpr(expname,station,startdate,enddate,vobs):
def readExpr(expname,station,startdate,enddate,vobsm):
for model in expname:
dmod = pd.read_csv(expname[model]["path"])
date = map(lambda x : str(x), dmod["datetime"].values)
vm = pd.Series(dmod[station].values,index=date)
vmod = vm.loc[startdate:enddate].values
vmod = vmod + biasMean(vobs,vmod)
lrms = np.sqrt(((vmod-np.mean(vmod)-vobs+ np.mean(vobs))**2).mean())
lcorr = np.corrcoef(vmod, vobs)[0,1]
modstd = np.std(vmod)
obstd = np.std(vobs)
vmod = vmod + biasMean(vobsm,vmod)
lcorr = np.ma.corrcoef(vmod, vobsm)[0,1]
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)
rms.append(round(lrms /obstd,3))
cor.append(round(lcorr,3))
print "sation", station," GEBCO", "rms", round(lrms/obstd,3), " cor", round(lcorr,3)
clegend.append(model)
vssh.append(vmod)
cstyle.append(expname[model]["style"])
cstyle.append(expname[model]["line"])
ccolor.append(expname[model]["color"])
def readOper(oper,station,startdate,enddate,vobs):
def readOper(oper,station,startdate,enddate,vobsm):
param = {'from': startdate, 'too': enddate}
for model in oper:
res1 = req.get("http://oceandata.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
voper = voper + biasMean(vobs,voper)
lrms = np.sqrt(((voper-np.mean(voper)-vobs+ np.mean(vobs))**2).mean())
lcorr = np.corrcoef(voper, vobs)[0,1]
modstd = np.std(voper)
obstd = np.std(vobs)
voper = voper + biasMean(vobsm,voper)
lcorr = np.ma.corrcoef(voper, vobsm)[0,1]
lrms = np.sqrt(np.mean((voper - vobsm)**2),dtype=np.float64)
lrms = np.sqrt(((voper-np.mean(voper)-vobsm+ np.mean(vobsm))**2).mean())
obstd = np.std(vobsm)
rms.append(round(lrms /obstd,3))
cor.append(round(lcorr,3))
print "sation", station," OPER", "rms", round(lrms /obstd,3), " cor", round(lcorr,3)
clegend.append(oper[model]["title"])
vssh.append(voper)
cstyle.append(oper[model]["style"])
cstyle.append(oper[model]["line"])
ccolor.append(oper[model]["color"])
##############################
def main():
dovalidation, station, startdate, enddate, obstyle, expname, oper,tickint = readConfig()
vobs = readObs(station, startdate, enddate,obstyle)
if dovalidation["experiment"] :readExpr(expname,station,startdate,enddate,vobs)
if dovalidation["operational"] :readOper(oper,station,startdate,enddate,vobs)
vobsm = readObs(station, startdate, enddate,obstyle)
if dovalidation["experiment"] :readExpr(expname,station,startdate,enddate,vobsm)
if dovalidation["operational"] :readOper(oper,station,startdate,enddate,vobsm)
stick = getDateRange(startdate,enddate)
itick = range(len(stick))
......@@ -161,7 +170,7 @@ def main():
ax.set_ylabel ('SSH(m)',fontsize='20',weight='bold')
for im in range(len(vssh)):
ax.plot(vssh[im],cstyle[im],lw=0.75)
ax.plot(vssh[im],linestyle=cstyle[im],color=ccolor[im],lw=0.75)
ax.legend(clegend, numpoints=1, prop=dict(size='small'),loc="best")
ax.grid()
......@@ -173,5 +182,14 @@ def main():
plt.savefig("ssh"+station+'.png', bbox_inches='tight',dpi=300,facecolor='w',edgecolor='w',orientation='portrait')
plt.close(1)
print ccolor[1:]
print clegend[1:]
print cor[:]
print ccolor[:]
print clegend[1:]
TaylorDiagram(RMSVEC=rms[:], CORVEC=cor[:],COLORVEC=ccolor[1:],LABELVEC=clegend[1:], station=station, info=startdate+"-"+enddate)
if __name__ == "__main__":
main()
taylorDiagram.py
View file @ e216a845
import numpy as np
import matplotlib.pyplot as plt
import math
#matplotlib.rcParams['xtick.direction'] = 'out'
#matplotlib.rcParams['ytick.direction'] = 'out'
#plt.tick_params(axis="y", labelcolor="b")
X=np.arange(0.0,1.1,0.01)
Y=np.arange(0.0,1.1,0.01)
print X.shape
print Y.shape
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])
fig=plt.figure(num=1,figsize=(3.0,3.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('RMS/OBS_STD',fontsize='5',weight='bold')
ax.set_ylabel ('RMS/OBS_STD',fontsize='5',weight='bold')
ax.grid(False)
vc=np.arange(0.0,1.2,0.2)
ax.contour(X,Y,h,vc,colors='0.5',linewidths=0.2)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_xlim(0.0,1.005)
ax.set_ylim(0.0,1.005)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
radius = np.arange(0.0,1.20,0.20)
xangle = [1.0,0.99,0.95] + list(np.arange(0.9,-0.1,-0.1))
xangle = list(np.arange(0.0,1.10,0.10))
#print xangle
for ang in xangle:
print ang, math.degrees(math.acos(ang))
nline = 0
for rd in radius:
for ang in xangle:
#print math.cos(math.acos(ang))
if ang == 0.99 or ang == 0.95:
ax.plot([0.0,rd*math.cos(math.acos(ang))],[0.0,rd*math.sin(math.acos(ang))],'k-',lw="0.0")
else:
#ax.plot([0.0,rd*math.cos(math.acos(ang))],[0.0,rd*math.sin(math.acos(ang))],'k-',lw="0.25")
#ax.plot([0.0,0.0],[rd*math.cos(math.acos(ang)),rd*math.sin(math.acos(ang))],'k-',lw="0.25")
#ax.plot([0.0,1.0],[0.0,1.0],'k-',lw="0.25")
if nline == 0:
nline += 1
ax.plot([0.0,rd*ang],[0.0,rd*math.sqrt(1.0 - (ang * ang))],color="0.5",ls="-",label="viken",lw="0.2")
if rd == 1.0: ax.annotate(str(ang), xy=(rd*ang, rd*math.sqrt(1.0 - (ang * ang))))
else:
nline += 1
ax.plot([0.0,rd*ang],[0.0,rd*math.sqrt(1.0 - (ang * ang))],color="0.5",ls="-",lw="0.2")
if rd == 1.0 and ang > 0.0: ax.annotate(str(ang), xy=(rd*ang, rd*math.sqrt(1.0 - (ang * ang))))
ax.plot(0.2*math.cos(math.acos(0.5)),0.2*math.sin(math.acos(0.5)),'ob',ms=3)
#ax.legend(numpoints=1,prop=dict(size='xx-small'),loc = 'upper right', bbox_to_anchor = (0.5, 0.5))
ax.legend(numpoints=1,prop=dict(size='xx-small'),loc = 'best')
plt.savefig("taylor.png", bbox_inches='tight',dpi=300,facecolor='w',edgecolor='w',orientation='portrait')
plt.close(1)
def TaylorDiagram(RMSVEC, 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)
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])
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('RMS/OBS_STD',fontsize='15',weight='bold',color="green")
ax.set_ylabel('RMS/OBS_STD',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)
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.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]
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)
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)
lenm = len(RMSVEC)
print "lenm", lenm
for ii in range(lenm):
vrms, vcor, vcol, vlab = RMSVEC[ii], CORVEC[ii], COLORVEC[ii],LABELVEC[ii]
print vrms, vcor, vcol, vlab
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)
ax.set_title(station+" "+info,fontsize="20")
plt.savefig("taylor"+station+".pdf", bbox_inches='tight',dpi=300,facecolor='w',edgecolor='w',orientation='portrait')
plt.close(1)
def main():
rms = [0.4,2.4,2.1,2.9]
cor = [0.5,0.25,0.76,0.24]
lab = ["exp1","exp2","exp3","exp4"]
col = ["b","g","r","k"]
TaylorDiagram(RMSVEC=rms, CORVEC=cor,COLORVEC=col,LABELVEC=lab, station="viken")
if __name__ == "__main__":main()
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