sshval.py 3.32 KB
Newer Older
saeed's avatar
saeed committed
1
2
3
4
5
6
7
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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
import os,sys
import numpy as np 
import pandas as pd
import requests as req
import matplotlib.pyplot as plt
import json
import datetime

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 fill_obs(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"] = -999.0
            
        initd = initd + step
    return dmiss
    
def bias_mean(vobs, vmodel):
    print "model shape", vmodel.shape[0]
    print "obs shape", vobs.shape[0]
    if vobs.shape[0] != vmodel.shape[0] :
        print "vobs and vmodel should have the same size"
        sys.exit()
    return np.mean(vobs- vmodel)

##############################
station = "barseback"
sfrom   = '2017100100'
stoo    = '2017102900'
param   = {'from': sfrom, 'too': stoo}

#read observation
res    = req.get("http://oceandata.smhi.se/ssh/"+station+"/OBSERVATION",params = param)
#res   = req.get("http://oceandata.smhi.se/ssh/"+station+"/OBSERVATION?"+"from="+sfrom+"&too="+stoo)
#print res.content
dobs     = res.json()
print len(dobs.keys())
date_obs = dobs.keys()
dmiss=fill_obs(date_obs,sfrom, stoo)
if dmiss.keys():
    dobs.update(dmiss)
print len(dobs.keys())
vo      = pd.DataFrame.from_dict(dobs,orient="index")
vob     = vo.loc[sfrom:stoo].values
vobs   = np.ma.masked_where(np.array(vob) == -999.0, np.array(vob))

#read model data 1 experiment
dmod1  = pd.read_csv("ssh_ns6_minh_4.csv")
date   = map(lambda x : str(x), dmod1["datetime"].values)
v1     = pd.Series(dmod1[station].values,index=date)
vmod1  = v1.loc[sfrom:stoo].values
vmod1  = vmod1 + bias_mean(vobs,vmod1)

#read model data 2 experiment
dmod2  = pd.read_csv("ssh_ns6_minh_4_correct_sss.csv")
date   = map(lambda x : str(x), dmod2["datetime"].values)
v2     = pd.Series(dmod2[station].values,index=date)
vmod2  = v2.loc[sfrom:stoo].values
vmod2  = vmod2 + bias_mean(vobs,vmod2)

##read operational
res1   = req.get("http://oceandata.smhi.se/ssh/"+station+"/NS01_NEMO",params = param)
doper  = res1.json()
vm     = pd.DataFrame.from_dict(doper,orient="index")
voper  = vm["raw"].loc[sfrom:stoo].values
voper  = voper + bias_mean(vobs,voper)

###plotting
ldate  = dict(v2.loc[sfrom:stoo]).keys()
idate  = [int(vd) for vd in ldate]
sidate = sorted(idate)

itick = range(len(sidate))
ctick = sidate
stick = 48

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')
ax.plot(vobs,'r-',vmod2,'b-',voper,'g-')
ax.legend(["OBS",'NEMO_NS01_GEBCO_northsea_6m_mindepth_4m',"NEMO_NS01_OPER"], loc="best")
#ax.set_ylim(-1.0,1.0)
ax.grid()
ax.set_title(station,fontsize="30")
ax.xaxis.set_ticks(itick[::stick])
ax.xaxis.set_ticklabels(ctick[::stick],rotation="vertical",fontsize="15")
plt.savefig("ssh"+station+'.png', bbox_inches='tight',dpi=300,facecolor='w',edgecolor='w',orientation='portrait')
plt.close(1)