-
Notifications
You must be signed in to change notification settings - Fork 2
/
Reanalysis_data.py
168 lines (118 loc) · 11.4 KB
/
Reanalysis_data.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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jul 17 12:03:32 2020
@author: anjakatzenberger
"""
import netCDF4 as nc
import subprocess
import os
import matplotlib.pyplot as plt
import statistics as stat
import numpy as np
from pyts.decomposition import SingularSpectrumAnalysis
###############################################################
# REANALYSIS DATA #
# anja.katzenberger@pik-potsdam.de #
# 18.06.2020 #
###############################################################
# Extract reanalysis data for pr for years of interest
# merges all pr files to one timeseries
# calculate monthly mean
# extract JJAS
# cutting India
# calculates averages over all grid cells
# Directories with files of interest
dir = ""
# Where to save
outdir = ""
mask = "fx_sftlf_CNRM-CM6-1-HR_masked_lonlat_missval.nc "
#files = os.listdir(dir)
# creates list of files of pr reanalysis data
#pr_reanalysis = []
#pr_reanalysis_dir = []
#periods = []
#for file in files:
# dir_file = os.sep.join((dir, file))
# data = nc.Dataset(dir_file)
# if file.startswith("gswp3-w5e5_obsclim_pr_global_daily"):
# period = file.split('_daily_')[1][:9]
# pr_reanalysis.append(file)
# periods.append(period)
# pr_reanalysis_dir.append(dir_file)
# landmask
mul_output = os.sep.join((outdir, "reanalysis_missmask.nc"))
cdo_cmd_mul = ' '.join(("cdo mul", "reanalysis.nc", mask, mul_output))
# calculates monthly mean
monmean_output = os.sep.join((outdir, "reanalysis_missmask_monmean.nc"))
cdo_cmd_monmean = ' '.join(("cdo -O monmean", mul_output, monmean_output))
# calculate JJAS mean
JJASmean_output = os.sep.join((outdir, "reanalysis_missmask_monmean_JJAS.nc"))
cdo_cmd_JJASmean = "cdo -O timselmean,4,5,8 " + monmean_output + " " + JJASmean_output #
# cut India
India_output = os.sep.join((outdir, "reanalysis_missmask_monmean_JJAS_India.nc"))
nco_cmd_India = "ncks -Od lon,67.5,98.0 -d lat,6.0,36.0 " + JJASmean_output + " " + India_output
# Calculate mean over India
output_fldmean = os.sep.join((outdir, "reanalysis_mask_monmean_JJAS_India_fldmean.nc"))
cdo_cmd_fldmean = "cdo -O fldmean " + India_output + ' ' + output_fldmean
#masking
#output_mask = os.sep.join((outdir, "reanalysis_monmean_JJAS_India_fldmean_mask.nc"))
#cdo_cmd_mask = "cdo -O mul landseamask.nc " + output_fldmean + ' ' + output_mask
#subprocess.check_call(cdo_cmd_mul, shell=True)
#subprocess.check_call(cdo_cmd_monmean, shell=True)
#subprocess.check_call(cdo_cmd_JJASmean, shell=True)
#subprocess.check_call(nco_cmd_India, shell=True)
#subprocess.check_call(cdo_cmd_fldmean, shell=True)
# extract pr data
data = nc.Dataset("reanalysis_monmean_JJAS_mul_India_fldmean.nc")
pr_data = []
for i in range(116):
pr = data['pr'][i,0,0] # [time, lat, lon]
pr = pr * 86400
pr_data.append(pr)
print(pr_data)
#----------------------------
### Results
#----------------------------
# results- without mask
#pr_data = [4.910306737292558, 6.455029989592731, 6.440762942656875, 5.890526412986219, 6.039390270598233, 6.307325349189341, 5.752819078043103, 6.230707000941038, 6.1676968820393085, 6.991062965244055, 6.4960206393152475, 6.0876826057210565, 6.287586595863104, 7.0048145251348615, 5.987271387130022, 6.475490424782038, 6.967211631126702, 6.13888178486377, 6.22219517827034, 5.6447421899065375, 6.839395361021161, 6.725191720761359, 6.745799258351326, 7.059773616492748, 5.849996558390558, 6.629764381796122, 6.966918054968119, 6.7373961908742785, 6.883983104489744, 6.762229464948177, 6.747045228257775, 6.675651529803872, 7.015927671454847, 6.941897445358336, 6.9533884059637785, 7.176703680306673, 6.737322639673948, 7.438399479724467, 6.912340549752116, 6.551996874623001, 6.770780892111361, 6.843538745306432, 6.593120796605945, 6.026717461645603, 7.010741368867457, 6.553406291641295, 6.592639884911478, 6.648014509119093, 6.83770808391273, 6.671668449416757, 5.829026293940842, 6.288261758163571, 6.534526264294982, 6.558050075545907, 6.470407219603658, 6.71029225923121, 6.302806036546826, 7.117875292897224, 7.4299115454778075, 6.6357823787257075, 7.582089607603848, 6.763324560597539, 6.804410135373473, 7.144559919834137, 6.090275757014751, 6.217145919799805, 6.723236641846597, 6.4443656941875815, 6.878781085833907, 7.083611120469868, 6.601673481054604, 5.658786068670452, 7.356878346763551, 6.7708255257457495, 7.236818899400532, 6.851535709574819, 6.836819811724126, 7.21993544138968, 6.0185488779097795, 6.975998170673847, 6.937028607353568, 6.578088691458106, 7.249260996468365, 6.747580203227699, 6.642046174965799, 6.284980243071914, 6.330348132178187, 7.622387493029237, 6.494440860114992, 6.626441376283765, 6.639065151102841, 6.6050587221980095, 6.664836988784373, 7.247607037425041, 7.021061168052256, 7.3055025190114975, 7.054129033349454, 7.2155594592913985, 6.603626045398414, 6.625077221542597, 6.355777359567583, 5.926699144765735, 6.854561367072165, 6.646568002179265, 7.071985001675785, 6.934010493569076, 7.639988232403994, 6.9651949452236295, 6.310125323943794, 7.609228114597499, 7.601571246050298, 6.855403748340905, 7.130786357447505, 6.393951061181724, 6.753709469921887, 6.754203583113849]
#plt_reanalysis = plt.plot(np.arange(1900, 2016, 1), pr_data_mask)
#mean_1900_1920 = stat.mean(pr_data[0:20]) # 6.224876229534857
#std_1900_1920 = stat.stdev(pr_data[0:20]) # 0.48735094106852334
#mean_1995_2015 = stat.mean(pr_data[95:115]) # 6.89267270732671
#std_1995_2015 = stat.stdev(pr_data[95:115]) # 0.4598834165650233
#mean_1900_2015 = stat.mean(pr_data[0:115]) # 6.6789324189860215
#std_1900_2015 = stat.stdev(pr_data[0:115]) # 0.46686049771979
# results - with mask (but 0_1)
#pr_data_mask = [3.4773014718666673, 4.2265024851076305, 4.103496274910867, 3.542446461506188, 3.8898721570149064, 4.089379473589361, 3.686978027690202, 4.027499654330313, 4.208183835726231, 4.512303590308875, 4.071065224707127, 4.022319638170302, 3.859594836831093, 4.3406891520135105, 3.946462890598923, 4.623491945676506, 4.565412586089224, 3.9602858014404774, 4.205869173165411, 3.5405419883318245, 4.4869777746498585, 4.229091864544898, 3.9572299690917134, 4.41015197429806, 3.764609433710575, 4.058751056436449, 4.018848587293178, 3.7981949863024056, 3.760857693850994, 3.7358125671744347, 4.311617254279554, 4.140651260968298, 4.212394484784454, 4.066775995306671, 4.042239126283675, 4.210812505334616, 3.948727890383452, 4.5153418206609786, 3.8966970168985426, 3.818390448577702, 3.6336474353447556, 4.369600431527942, 3.980747179593891, 4.074188007507473, 4.1289204731583595, 4.034238704480231, 4.289600299671292, 4.277937090955675, 4.203458642587066, 3.9891584194265306, 3.4927943721413612, 3.86222853558138, 4.244584450498223, 4.207520303316414, 4.180410399567336, 4.363775427918881, 3.7339398404583335, 4.186864674557, 4.512681404594332, 4.028070461936295, 4.6478021889925, 3.978953661862761, 4.033759678713977, 4.430002940353006, 3.6127111176028848, 3.6543684429489076, 3.9468284463509917, 3.739823622163385, 4.045426030643284, 4.348901112098247, 4.167301312554628, 3.3113803365267813, 4.320868360809982, 3.955947223585099, 4.42898862529546, 4.088593041524291, 4.058859497308731, 4.324464511591941, 3.554136073216796, 4.2196483933366835, 4.049274267163128, 3.7939485046081245, 4.245328134857118, 4.102368489839137, 4.006647574715316, 3.663677698932588, 3.7650570273399353, 4.478833707980812, 4.024416476022452, 4.302540596108884, 4.042615368962288, 3.7470913608558476, 4.085420595947653, 4.323129274416715, 4.146021127235144, 4.079543729312718, 4.167315142694861, 4.113991151098162, 4.015892080496997, 3.906021674629301, 3.815943142399192, 3.591736138332635, 4.087749402970076, 3.9101210539229214, 4.053009976632893, 3.933431440964341, 4.45376563584432, 4.267013794742525, 3.5849285661242902, 4.384613363072276, 4.409501957707107, 4.063673014752567, 4.22739798668772, 3.68373611709103, 3.897093376144767, 3.985403850674629]
#mean_1900_1920 = stat.mean(pr_data_mask[0:20]) # 4.044984833453782
#std_1900_1920 = stat.stdev(pr_data_mask[0:20]) # 0.3273863320688452
#mean_1995_2015 = stat.mean(pr_data_mask[95:115]) # 4.032323937281035
#std_1995_2015 = stat.stdev(pr_data_mask[95:115]) # 0.25032180669647686
#mean_1900_2015 = stat.mean(pr_data_mask[0:115]) # 4.055390004372305
#std_1900_2015 = stat.stdev(pr_data_mask[0:115]) # 0.2744362597459318
# results with missing values mask
pr_data_mask = [5.251554446294904, 6.383026507683098, 6.197258178144693, 5.34993892069906, 5.87463432457298, 6.175937759689987, 5.568216252140701, 6.082484987564385, 6.355361198075116, 6.8146544974297285, 6.148279365152121, 6.074662157334387, 5.828908737748861, 6.555475783534348, 5.960100190714002, 6.9825750309973955, 6.894861767068505, 5.9809761587530375, 6.351865315809846, 5.3470628801733255, 6.776406615972519, 6.386937294155359, 5.976361292414367, 6.6603811690583825, 5.685458122752607, 6.129682227037847, 6.069419905543327, 5.7361801620572805, 5.679792165756226, 5.641967989504337, 6.511570117436349, 6.253370828926563, 6.3617205480113626, 6.141801201738417, 6.104744598269463, 6.359331076964736, 5.963520635850728, 6.819242960773408, 5.88494217954576, 5.7666800217702985, 5.48767454456538, 6.599138793535531, 6.011877721175551, 6.152995442971587, 6.235654419288039, 6.0926627134904265, 6.478319317102432, 6.46070537623018, 6.348224845714867, 6.0245807049795985, 5.274952528998256, 5.832886160351336, 6.410334748215973, 6.354359141550958, 6.313416897319257, 6.5903415670618415, 5.639139725826681, 6.323164002969861, 6.8152246763929725, 6.083347485400736, 7.019289652816951, 6.009168899618089, 6.091939145699143, 6.690361141227186, 5.456055072136223, 5.518967751413584, 5.9606521390378475, 5.648025590926409, 6.1095574870705605, 6.567877647466958, 6.293619051575661, 5.000973679125309, 6.525541702285409, 5.974423815496266, 6.688829138875008, 6.174750882200897, 6.129845674149692, 6.530972546897829, 5.367593094706535, 6.372675276361406, 6.115369917824864, 5.729766748845577, 6.411458132788539, 6.195555184967816, 6.050993129611015, 5.5330273462459445, 5.686133913695812, 6.764107220806181, 6.077828630805016, 6.497861933894455, 6.105313519947231, 5.6590016931295395, 6.169959367252886, 6.528956489637494, 6.2614803202450275, 6.161083560436964, 6.29363979678601, 6.213108147494495, 6.064954656176269, 5.899024405516684, 5.7629842311143875, 5.424378393217921, 6.1734766233712435, 5.905215279199183, 6.121011357754469, 5.940419901162386, 6.726247840560973, 6.444207904860377, 5.414096941240132, 6.621811422519386, 6.6593998577445745, 6.137115298770368, 6.384379346854985, 5.563320382498205, 5.88554001878947, 6.018910347484052]
plt_reanalysis = plt.plot(np.arange(1900, 2016, 1), pr_data_mask,marker=".")
stat.mean(pr_data_mask) # 6.1
# values
mean_1900_1930 = stat.mean(pr_data_mask[0:30]) # 6.097347380127758
std_1900_1930 = stat.stdev(pr_data_mask[0:30]) # 0.46266551278823365
mean_1985_2015 = stat.mean(pr_data_mask[85:115]) # 6.102636193390936
std_1985_2015 = stat.stdev(pr_data_mask[85:115]) # 0.38143934169677174
mean_1900_2015 = stat.mean(pr_data_mask[0:115]) # 6.124606040792297
std_1900_2015 = stat.stdev(pr_data_mask[0:115]) # 0.4144641707853791
# singular spectrum analysis
L = 20
F = pr_data_mask
F_arr = np.array(F)
F_in = np.array([F_arr])
# Singular Spectrum Analysis
ssa = SingularSpectrumAnalysis(window_size=L)
X_ssa = ssa.transform(F_in)
data_ssa = X_ssa[0]
# 20 years averaging
period = 20
timeseries_av = []
for m in range(period, len(pr_data_mask)):
timeseries_av.append(stat.mean(pr_data_mask[m : m + period]))