-
Notifications
You must be signed in to change notification settings - Fork 0
/
obsVSmodel_deepestbottom.py
122 lines (119 loc) · 4.76 KB
/
obsVSmodel_deepestbottom.py
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
'''
If the deepest observation depth is >50m(or <50m, or all), draw the correlation of this observation and appriate model data.
'''
import numpy as np
import pandas as pd
from datetime import datetime, timedelta
import netCDF4
import matplotlib.pyplot as plt
# import sys
# sys.path.append('../moj')
# import utilities
# from matplotlib import path
from scipy import stats
from watertempModule import water_roms
from turtleModule import mon_alpha2num, np_datetime, mean_value, bottom_value, index_by_depth
def show2pic(x1, y1, fontsize):
FONTSIZE = fontsize
fig1 = plt.figure()
ax1 = fig1.add_subplot(111)
x = np.arange(0.0, 30.0, 0.01)
'''
for i in range(10):
# ax.plot(temp[index[i]], obstemp[index[i]], '.', color=colors[i], label='{0}'.format(i))
ax.scatter(temp[index[i]], obsdata['temp'][index[i]], s=50, c=colors[i], label='{0}'.format(i))
'''
# ax.scatter(temp[index[0]], obsdata['temp'][index[0]], s=50, c='b', label='<45')
# ax.scatter(temp[index[1]], obsdata['temp'][index[1]], s=50, c='r', label='>=45')
ax1.scatter(x1, y1, s=50, c='b')
ax1.plot(x, x, 'r-', linewidth=2)
plt.axis([0, 30, 0, 30], fontsize=15)
plt.xlabel('Model temp', fontsize=FONTSIZE)
plt.ylabel('CTD temp', fontsize=FONTSIZE)
i = x1[x1.isnull()==False].index
fit = np.polyfit(x1[i], y1[i], 1)
fit_fn = np.poly1d(fit)
x2, y2 = x1[i], fit_fn(x1[i])
plt.plot(x2, y2,'y-', linewidth=2)
gradient, intercept, r_value, p_value, std_err = stats.linregress(y1[i], x1[i])
r_squared = r_value**2
# ax1.set_title('R-squard: %.4f' % r_squared, fontsize=FONTSIZE)
plt.savefig('obsVSmodelDeepestBottom1.png',dpi=200)
fig2 = plt.figure()
ax2 = fig2.add_subplot(111)
nbins = 200
H, xedges, yedges = np.histogram2d(x1[i], y1[i], bins=nbins)
H = np.rot90(H)
H = np.flipud(H)
Hmasked = np.ma.masked_where(H==0, H)
plt.pcolormesh(xedges, yedges, Hmasked)
plt.xlabel('Model temp', fontsize=FONTSIZE)
plt.ylabel('CTD temp', fontsize=FONTSIZE)
cbar = plt.colorbar()
cbar.ax.set_ylabel('Counts', fontsize=FONTSIZE)
cbar.ax.set_yticks(fontsize=20)
plt.axis([0, 30, 0, 30])
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.plot(x, x, 'r-', linewidth=2)
plt.plot(x2, y2, 'y-', linewidth=2)
# plt.title('R-squard: %.4f' % r_squared, fontsize=FONTSIZE)
plt.savefig('obsVSmodelDeepestBottom2.png',dpi=200)
return ax1, ax2, r_squared
#############################MAIN CODE###########################################
FONTSIZE = 25
# obs = pd.read_csv('ctd_extract_TF.csv')
obs = pd.read_csv('ctd_good.csv')
tf_index = np.where(obs['TF'].notnull())[0]
obslat, obslon = obs['LAT'][tf_index].values, obs['LON'][tf_index].values
obstime = np_datetime(obs['END_DATE'][tf_index])
# obsdepth = mean_value(obs['TEMP_DBAR'][tf_index])
obsdepth = obs['MAX_DBAR'][tf_index].values
# obstemp = mean_value(obs['TEMP_VALS'][tf_index])
obstemp = bottom_value(obs['TEMP_VALS'][tf_index])
obsdata = pd.DataFrame({'depth':obsdepth, 'temp':obstemp, 'lon':obslon,
'lat':obslat, 'time':obstime}).sort_index(by='depth')
starttime = datetime(2009, 8, 24)
endtime = datetime(2013,12 ,13)
tempobj = wtm.water_roms()
url = tempobj.get_url(starttime, endtime)
# temp = tempobj.watertemp(obslon, obslat, obsdepth, obstime, url)
temp = tempobj.watertemp(obsdata['lon'].values, obsdata['lat'].values,
obsdata['depth'].values, obsdata['time'].values, url)
temp = pd.Series(temp, index = obsdata['temp'].index)
index = index_by_depth(obsdata['depth'], 50)
# colors = utilities.uniquecolors(10)
tp='all'
if tp == 'all':
x1, y1 = temp, obsdata['temp']
ax1, ax2, r_squared = show2pic(x1, y1, FONTSIZE)
ax1.set_title('R-squared: %.4f' % r_squared, fontsize=FONTSIZE)
ax2.set_title('R-squared: %.4f' % r_squared, fontsize=FONTSIZE)
elif tp == '<50':
x1, y1 = temp[index[0]], obsdata['temp'][index[0]]
ax1, ax2, r_squared = show2pic(x1, y1, FONTSIZE)
ax1.set_title('%s, R-squared: %.4f' % (tp, r_squared), fontsize=FONTSIZE)
ax2.set_title('%s, R-squared: %.4f' % (tp, r_squared), fontsize=FONTSIZE)
elif tp == '>50':
x1, y1 = temp[index[1]], obsdata['temp'][index[1]]
ax1, ax2, r_squared = show2pic(x1, y1, FONTSIZE)
ax1.set_title('%s, R-squared: %.4f' % (tp, r_squared), fontsize=FONTSIZE)
ax2.set_title('%s, R-squared: %.4f' % (tp, r_squared), fontsize=FONTSIZE)
plt.show()
'''
# Plot Deepest Data Quantity
fig = plt.figure()
ax = fig.add_subplot(111)
y = obsdata['depth'].values
x = np.arange(1, np.amax(y)+1)
bar = np.array([0]*np.amax(y))
for i in y:
if i in x:
bar[i-1] = bar[i-1]+1
plt.barh(x, bar)
plt.ylim([250, 0])
plt.ylabel('depth', fontsize=25)
plt.xlabel('Quantity', fontsize=25)
plt.title('Deepest data histogram', fontsize=25)
plt.show()
'''