sshval.py 6.81 KB
Newer Older
saeed's avatar
saeed committed
1
2
3
4
5
6
7
import os,sys
import numpy as np 
import pandas as pd
import requests as req
import matplotlib.pyplot as plt
import json
import datetime
8
from taylorDiagram import TaylorDiagram
saeed's avatar
saeed committed
9
10
11
12
reload(sys)
sys.setdefaultencoding('utf8')

clegend = []
saeed's avatar
saeed committed
13
vssh    = []        
saeed's avatar
saeed committed
14
cstyle  = [] 
15
ccolor  = [] 
saeed's avatar
saeed committed
16
17
18
rms     = []
cor     = []  

19
20
21
22
23
24
25
26
27
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)

saeed's avatar
saeed committed
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
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])


def getDateRange(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)
    dateticks = []
    while initd <= totd:
        dateticks.append(strDate(initd)) 
        initd = initd + step
    
    return dateticks

44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

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

saeed's avatar
saeed committed
61
62
63
64
65
66
def biasMean(vobs, vmodel):
    
    if vobs.shape[0] != vmodel.shape[0] :
        print "vobs and vmodel should have the same size"
        sys.exit()
    
67
    return np.mean(vobs- vmodel)
saeed's avatar
saeed committed
68

saeed's avatar
saeed committed
69
70
71
72
73
74
75
76
77
78
79
80
81
82

def readConfig():

    if len(sys.argv) != 2:
        print "usage "+sys.argv[0]+" config file"
        sys.exit(1)
        
    config_file = sys.argv[1]
    if not os.path.isfile(config_file):
        print 'specified config file {} does not exist.'.format(config_file)
        sys.exit(1)

    with open(config_file, 'r') as fptr_cfg:
        cfg = json.load(fptr_cfg)
83
84
        dovalidation = cfg["dovalidation"]
        station      = cfg["station"]
85
86
        startdate    = format_time_string(cfg["startdate"])
        enddate      = format_time_string(cfg["endate"])
87
88
89
90
91
        obstyle      = cfg["obstyle"]
        expname      = cfg["experiment"]
        oper         = cfg["operational"]
        tickint      = cfg["tickinterval_hour"]
        return dovalidation, station, startdate, enddate,obstyle, expname, oper, tickint 
saeed's avatar
saeed committed
92
93
    

saeed's avatar
saeed committed
94
95
96
def readObs(station, startdate,enddate,obstyle):
    param   = {'from': startdate, 'too': enddate}
    res    = req.get("http://oceandata.smhi.se/ssh/"+station+"/OBSERVATION",params = param)
97
    dobs   = res.json()
saeed's avatar
saeed committed
98
    date_obs = dobs.keys()
99
    dmiss = fillObs(date_obs,startdate, enddate)
saeed's avatar
saeed committed
100
101
    if dmiss.keys():
        dobs.update(dmiss)
102
103
104
105
106
    
    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)
saeed's avatar
saeed committed
107
    clegend.append("OBSERVATION")
108
109
110
    vssh.append(vobsm)
    cstyle.append(obstyle["line"]) 
    ccolor.append(obstyle["color"]) 
saeed's avatar
saeed committed
111
        
112
    return vobsm
saeed's avatar
saeed committed
113
    
114
def readExpr(expname,station,startdate,enddate,vobsm):
saeed's avatar
saeed committed
115
        
saeed's avatar
saeed committed
116
    for model in expname:
saeed's avatar
saeed committed
117
118
        dmod   = pd.read_csv(expname[model]["path"])
        date   = map(lambda x : str(x), dmod["datetime"].values)
saeed's avatar
saeed committed
119
        vm     = pd.Series(dmod[station].values,index=date)
saeed's avatar
saeed committed
120
        vmod   = vm.loc[startdate:enddate].values
121
122
123
124
125
        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)
saeed's avatar
saeed committed
126
127
        rms.append(round(lrms /obstd,3))
        cor.append(round(lcorr,3))
128
        print "sation", station," GEBCO", "rms", round(lrms/obstd,3), " cor", round(lcorr,3)
saeed's avatar
saeed committed
129
        clegend.append(model)
130
        vssh.append(vmod)
131
132
        cstyle.append(expname[model]["line"]) 
        ccolor.append(expname[model]["color"]) 
saeed's avatar
saeed committed
133
        
134
def readOper(oper,station,startdate,enddate,vobsm):
saeed's avatar
saeed committed
135
136
137
138
139
140
    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
141
142
143
144
145
        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)
saeed's avatar
saeed committed
146
147
        rms.append(round(lrms /obstd,3))
        cor.append(round(lcorr,3))
148
        print "sation", station," OPER", "rms", round(lrms /obstd,3), " cor", round(lcorr,3)
saeed's avatar
saeed committed
149
        clegend.append(oper[model]["title"])
150
        vssh.append(voper)
151
152
        cstyle.append(oper[model]["line"])
        ccolor.append(oper[model]["color"])
saeed's avatar
saeed committed
153
154
##############################

saeed's avatar
saeed committed
155
156
def main():

157
    dovalidation, station, startdate, enddate, obstyle, expname, oper,tickint = readConfig()
saeed's avatar
saeed committed
158

159
160
161
    vobsm = readObs(station, startdate, enddate,obstyle)
    if dovalidation["experiment"]  :readExpr(expname,station,startdate,enddate,vobsm)
    if dovalidation["operational"] :readOper(oper,station,startdate,enddate,vobsm)
saeed's avatar
saeed committed
162
163
164

    stick = getDateRange(startdate,enddate)
    itick = range(len(stick))
saeed's avatar
saeed committed
165
166
    iitick=itick[::tickint]

saeed's avatar
saeed committed
167
168
169
170
171
    fig=plt.figure(num=1,figsize=(12.0,6.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('Time',fontsize='20',weight='bold')
    ax.set_ylabel ('SSH(m)',fontsize='20',weight='bold')

172
    for im in range(len(vssh)):
173
        ax.plot(vssh[im],linestyle=cstyle[im],color=ccolor[im],lw=0.75)
saeed's avatar
saeed committed
174

saeed's avatar
saeed committed
175
    ax.legend(clegend, numpoints=1, prop=dict(size='small'),loc="best")
saeed's avatar
saeed committed
176
177
178
    ax.grid()
    ax.set_title(station,fontsize="30")
    ax.xaxis.set_ticks(itick[::tickint])
179
    ax.set_xlim(iitick[0],iitick[-1])
saeed's avatar
saeed committed
180
    ax.xaxis.set_ticklabels(stick[::tickint],rotation="vertical",fontsize="15")
181
    print "Creating figure for station", station, " ssh"+station+'.png'
saeed's avatar
saeed committed
182
183
184
    plt.savefig("ssh"+station+'.png', bbox_inches='tight',dpi=300,facecolor='w',edgecolor='w',orientation='portrait')
    plt.close(1)

185
186
187
188
189
190
191
192
193

    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)

        
saeed's avatar
saeed committed
194
195
if __name__ == "__main__":
    main()