sshval.py 4.86 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
saeed's avatar
saeed committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
reload(sys)
sys.setdefaultencoding('utf8')


clegend = []
vexpr   = []        
cstyle  = [] 

def readConfig():

    if len(sys.argv) != 2:
        print "usage "+sys.argv[0]+" config file"
        sys.exit(1)
        
    config_file = sys.argv[1]
    """Read and parse config file"""
    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)
        station   = cfg["station"]
        startdate = cfg["startdate"]
        enddate   = cfg["endate"]
        obstyle   = cfg["obstyle"]
        expname   = cfg["experiment"]
        oper      = cfg["operational"]
        tickint   = cfg["tickinterval"]
        return station, startdate, enddate,obstyle, expname, oper, tickint 
saeed's avatar
saeed committed
38
39
40
41
42

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])

saeed's avatar
saeed committed
43
44
45
46
47
48
49
50
51
52
53
54
55

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
    
def fillObs(date_obs,fromt, tot):
saeed's avatar
saeed committed
56
57
58
59
60
61
62
63
64
65
66
67
    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"] = -999.0
            
        initd = initd + step
    return dmiss
    
saeed's avatar
saeed committed
68
def biasMean(vobs, vmodel):
saeed's avatar
saeed committed
69
    
saeed's avatar
saeed committed
70
71
72
73
74
    if vobs.shape[0] != vmodel.shape[0] :
        print "vobs and vmodel should have the same size"
        sys.exit()
    return np.mean(vobs- vmodel)

saeed's avatar
saeed committed
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
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()
    date_obs = dobs.keys()
    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.ma.masked_where(np.array(vob) == -999.0, np.array(vob))
    clegend.append("OBSERVATION")
    vexpr.append(vobs)
    cstyle.append(obstyle) 
        
    return vobs

def readExpr(expname,station,startdate,enddate,vobs):
saeed's avatar
saeed committed
93
        
saeed's avatar
saeed committed
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
    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)
        clegend.append(model)
        vexpr.append(vmod)
        cstyle.append(expname[model]["style"]) 
        
def readOper(oper,station,startdate,enddate,vobs):
    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)
        clegend.append(oper[model]["title"])
        vexpr.append(voper)
        cstyle.append(oper[model]["style"]) 
saeed's avatar
saeed committed
115
116
##############################

saeed's avatar
saeed committed
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def main():

    station, startdate, enddate, obstyle, expname, oper,tickint = readConfig()

    vobs = readObs(station, startdate, enddate,obstyle)
    readExpr(expname,station,startdate,enddate,vobs)
    readOper(oper,station,startdate,enddate,vobs)

    stick = getDateRange(startdate,enddate)
    itick = range(len(stick))
    
    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')

    for im in range(len(vexpr)):
        ax.plot(vexpr[im],cstyle[im])

    ax.legend(clegend, loc="best")
    ax.grid()
    ax.set_title(station,fontsize="30")
    ax.xaxis.set_ticks(itick[::tickint])
    ax.xaxis.set_ticklabels(stick[::tickint],rotation="vertical",fontsize="15")
saeed's avatar
saeed committed
141
    print "Creating figure ", "ssh"+station+'.png'
saeed's avatar
saeed committed
142
143
144
145
146
    plt.savefig("ssh"+station+'.png', bbox_inches='tight',dpi=300,facecolor='w',edgecolor='w',orientation='portrait')
    plt.close(1)

if __name__ == "__main__":
    main()