-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpp_run_initial_analysis
121 lines (108 loc) · 6.04 KB
/
pp_run_initial_analysis
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
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 22 10:02:45 2020
@author: josmarti
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import rms_plots as rpl
import nc as nc
import numpy as np
import re
#import cdo; c=cdo.Cdo()
metric='SIV'
syear=1979; eyear=2016;
obsin=np.zeros((12,40))
if metric=='SIE':
for m in ["%02d" % i for i in range(1,13)]:
obsin[int(m)-1,:]=np.loadtxt('/home/josmarti/Data/Observations/NSIDC/north/N_%s_extent_v3.0.csv' % m,usecols=4, skiprows=1, delimiter=',', max_rows=40)
obsin[obsin == -9999] = np.nan
if metric=='SIA':
for m in ["%02d" % i for i in range(1,13)]:
obsin[int(m)-1,:]=np.loadtxt('/home/josmarti/Data/Observations/NSIDC/north/N_%s_extent_v3.0.csv' % m,usecols=5, skiprows=1, delimiter=',', max_rows=40)
obsin[obsin == -9999] = np.nan
mean_obs=np.nanmean(obsin[:,(syear-1979):(eyear-1979)], axis=1)*1e12
annual_mean_obs=np.nanmean(obsin[:,(syear-1979):(eyear-1979)], axis=0)*1e12
gridpoint=nc.getvar('/home/josmarti/Data/areacella_fx_CanCM4_decadal2001_r0i0p0.nc','areacella')
PIOMAS=np.loadtxt('/home/josmarti/Data/Observations/PIOMAS_SIV_monthly_1979-2020.txt',usecols=range(1,13))
if metric=='SIV':
mean_obs=np.mean(PIOMAS[(syear-1979):(eyear-1979)], axis=0)*1e12
years=[2312, 2435, 2413, 2356, 2387, 2332]
ensembles=["%03d" % i for i in range(1,13)]
#ensembles=["001"]
run_length=1 #years
daily_matrix=np.zeros((6,12,1095))
days=365*3
monthly_matrix=np.zeros((6,12,36))
for e in ensembles:
for year in years:
if metric=='SIE':
obs=annual_mean_obs
geosum1=(nc.getvar('/home/josmarti/Data/PPrun/SIE/SIE_PPrun_e%s_i%s_%s.nc' % (e,str(year),str(year)), 'sicn')*np.sum(gridpoint[32:64])).squeeze()
geosum2=(nc.getvar('/home/josmarti/Data/PPrun/SIE/SIE_PPrun_e%s_i%s_%s.nc' % (e,str(year),str(year+1)), 'sicn')*np.sum(gridpoint[32:64])).squeeze()
geosum3=(nc.getvar('/home/josmarti/Data/PPrun/SIE/SIE_PPrun_e%s_i%s_%s.nc' % (e,str(year),str(year+2)), 'sicn')*np.sum(gridpoint[32:64])).squeeze()
if metric=='SIA':
obs=annual_mean_obs
var1=np.multiply(nc.getvar(('/home/josmarti/Data/PPrun/sc_dhfp1e_e%s_i%s_m01_%s_m01_%s_m12_sicn.nc' % (e,str(year),str(year),str(year))), 'sicn'), gridpoint)
var2=np.multiply(nc.getvar(('/home/josmarti/Data/PPrun/sc_dhfp1e_e%s_i%s_m01_%s_m01_%s_m12_sicn.nc' % (e,str(year),str(year+1),str(year+1))), 'sicn'), gridpoint)
var3=np.multiply(nc.getvar(('/home/josmarti/Data/PPrun/sc_dhfp1e_e%s_i%s_m01_%s_m01_%s_m12_sicn.nc' % (e,str(year),str(year+2),str(year+2))), 'sicn'), gridpoint)
if metric=='SIV':
obs=np.mean(PIOMAS[(syear-1979):(eyear-1979)], axis=1)*1e12
var1=np.multiply(nc.getvar(('/home/josmarti/Data/PPrun/sc_dhfp1e_e%s_i%s_m01_%s_m01_%s_m12_sic.nc' % (e,str(year),str(year),str(year))), 'sic'), gridpoint/913)
var2=np.multiply(nc.getvar(('/home/josmarti/Data/PPrun/sc_dhfp1e_e%s_i%s_m01_%s_m01_%s_m12_sic.nc' % (e,str(year),str(year+1),str(year+1))), 'sic'), gridpoint/913)
var3=np.multiply(nc.getvar(('/home/josmarti/Data/PPrun/sc_dhfp1e_e%s_i%s_m01_%s_m01_%s_m12_sic.nc' % (e,str(year),str(year+2),str(year+2))), 'sic'), gridpoint/913)
var=nc.getvar('/home/josmarti/Data/PPrun_Jan/sc_dhfp1e_e%s_i%s_m01_daily_sicn.nc4' % (e,str(year)),'sicn')
if metric != 'SIE':
var1nh=np.delete(var1, range(32), axis=1)
var2nh=np.delete(var2, range(32), axis=1)
var3nh=np.delete(var3, range(32), axis=1)
geosum1=np.sum(np.sum(var1nh, axis=1), axis=1)
geosum2=np.sum(np.sum(var2nh, axis=1), axis=1)
geosum3=np.sum(np.sum(var3nh, axis=1), axis=1)
daily_matrix[years.index(year),ensembles.index(e),:]=np.sum(np.sum(var[:,32:64,:], axis=1), axis=1)
if year==2413 and e=='011': #2415 e11 seems to have an extra timestep
geosum3=np.delete(geosum3, 12)
monthly_matrix[years.index(year),ensembles.index(e),:]=np.concatenate((geosum1,geosum2,geosum3))
#run=np.concatenate([geosum1, geosum2, geosum3])
run=geosum1
if len(run)>36:
print("%s, %s" % (str(year), e))
months=range(1,13)
plt.plot(months, geosum1, label=('%se%s' % (str(year),e)))
if run_length==2:
plt.plot(months, geosum2, label=('%se%s' % (str(year+1),e)))
if run_length==3:
plt.plot(months, geosum2, label=('%se%s' % (str(year+1),e)))
plt.plot(months, geosum3, label=('%se%s' % (str(year+2),e)))
#plt.plot(range(1,13), (np.mean(run/obsin[start:end])*obsin[start:end]), label='observations')
plt.plot(months, mean_obs, label='mean_obs', lw=2.5, color='black')
plt.title("Monthly mean %s for a %i year run of %i ensembles \n (observations: %i-%i)" % (metric,run_length,len(ensembles),syear,eyear))
if (len(years)*run_length*len(ensembles))<12:
plt.legend()
#plt.legend()
#%%
c=['g','r','c','m','y','k']
plt.figure(figsize=(10,7))
plt.plot(range(syear,eyear), obs)
if metric=='SIE':
old=nc.getvar('/home/josmarti/Data/Observations/NSIDC_1979_2010_nh_siea.nc', 'sie').squeeze()*1e12
plt.plot(range(syear,syear+32),np.mean(np.reshape(old,(12,32)), axis=0), label='old obs')
if metric=='SIA':
old=nc.getvar('/home/josmarti/Data/Observations/NSIDC_1979_2010_nh_siea.nc', 'sia').squeeze()*1e12
plt.plot(range(syear,syear+32),np.mean(np.reshape(old,(12,32)), axis=0), label='old obs')
plt.title('Observed annual mean %s with PPrun means' % metric)
for i in range(len(monthly_matrix)):
plt.hlines(np.mean(monthly_matrix, axis=(1,2))[i],syear,eyear, color=c[i], label=str(years[i]))
plt.legend()
#%%
"""
fig=plt.figure(figsize=(10,10))
for y in years:
fig.add_subplot(3,2,(years.index(y)+1), title="Std of %s" % str(y))
for i in range(12):
#plt.plot(range(days), daily_matrix[years.index(y),i,0:days], label='%s' % ensembles[i])
plt.plot(range(days), np.std(daily_matrix[years.index(y),:,0:days], axis=0))
#plt.title("First %i daily SIA values of %i for all 12 ensembles" % (days,y))
plt.legend()
fig.tight_layout()"""