-
Notifications
You must be signed in to change notification settings - Fork 0
/
cvdp_data.functions.ncl
2336 lines (2214 loc) · 89.3 KB
/
cvdp_data.functions.ncl
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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
;=================================================================================================
; These functions are not required to run the CVDP software.
; Rather they are aimed at operating on resultant cvdp_data.*.nc output files
; The aim of them is to allow easy plotting of say
; - change in climate mode between two runs
; - ensemble mean plots for certain experiments
; - difference with observations
;
; Author: Chris Brierley, started in Summer 2017
; Forked off and rejigged for Brown et al (2020) in Summer 2019
; Tidied up and moved to the past2future analyzer in Summer 2020
;
;
; NOTE: all programs assume filenames of the form exptname/modelname_exptname.cvdp_data.Y-Y.nc
; The cvdp_data files download from http://www.cesm.ucar.edu/working_groups/CVC/cvdp/data-repository.html
; will need to untarred in a directory of the experiment name and tidied up:
; > cd $data_dir/ssp585
; > wget http://webext.cgd.ucar.edu/Multi-Case/CVDP_repository/cmip6.ssp585/cmip6.ssp585.cvdp_data.tar
; > tar -xf cmip6.ssp585.cvdp_data.tar
; REMEMBER: these files use a subtly different version of CVDP.
; - Monsoon, IPCC regions and sea ice area diagnostics will be missing.
; - Less information about IOD, AMM, and ATL3
; - Definition of AMOC differs!
;=================================================================================================
; Check that the user-specified climatological period is within the time range of the data
; copied from functions.ncl
undef("check_custom_climo")
procedure check_custom_climo(mn:string,startyear:numeric,endyear:numeric,climo_startyear:numeric,climo_endyear:numeric)
local startyear,endyear,climo_startyear,climo_endyear,mn
begin
do gg = 0,dimsizes(startyear)-1
if (climo_startyear.ge.startyear(gg).and.climo_endyear.le.endyear(gg)) then
else
print("check_custom_climo: Warning! Beginning and/or ending of climatological period is outside time range of data.")
print("Dataset: "+mn+", years = "+startyear(gg)+":"+endyear(gg)+", set climatological period = "+climo_startyear+":"+climo_endyear)
print("The diagnostics package will proceed, but one or more dataset(s) will not have the full climatological period removed and/or tha package may fail with the following message: fatal:NclOneDValGetRangeIndex: start coordinate index out of range.")
end if
end do
end
;=================================================================================================
; find netcdf files that contain a desired variable and fit a filename string.
;
undef("find_files_wVar")
function find_files_wVar(name_str:string,variable_name:string)
local data_dir,name_str,variable_name,fil
begin
present_dir=systemfunc("pwd")
if ismissing(str_match_ind(present_dir,"scripts")) then
data_dir=present_dir+"/data/"
else
data_dir=present_dir+"/../data/"
end if
if isatt(variable_name,"ensemble") then
data_dir=data_dir+variable_name@ensemble+"/"
else
data_dir=data_dir+"{PMIP4,PMIP3}/"
end if
if isatt(name_str,"quiet") then
QUIET=name_str@quiet
else
QUIET=False
end if
ncfiles=systemfunc("ls "+data_dir+"/*.cvdp_data.[0-9]*-*.nc")
if any(ismissing(ncfiles)) then
if .not.QUIET then
print(variable_name)
print("find_files_wVar: There are no cvdp_data netcdf files in "+data_dir+". You may want to run the download_cvdp_data_data.sh script")
end if
ncfiles_wVar="missing"
else
if any(ismissing(str_match_ind(ncfiles,name_str))) then
if .not.QUIET then
print("find_files_wVar: There are no cvdp_data netcdf files in "+data_dir+"/ containing the string "+name_str)
end if
ncfiles_wVar="missing"
else
ncfiles_match=str_match(ncfiles,name_str)
relevant=new((/dimsizes(ncfiles_match),dimsizes(variable_name)/),logical)
relevant=(/False/)
do i=0,dimsizes(ncfiles_match)-1
this_file=ncfiles_match(i)
if isfilepresent(this_file) then
fil = addfile (this_file, "r")
do var_i=0,dimsizes(variable_name)-1
if isfilevar(fil,variable_name(var_i)) then
relevant(i,var_i) = True
end if
end do
end if
end do
if all(.not.(dim_num(relevant).eq.dimsizes(variable_name))) then
ncfiles_wVar="missing"
if .not.QUIET then
print("find_files_wVar: There are no cvdp_data netcdf files that contain "+name_str+" and have the requested variable")
end if
else
ncfiles_wVar=ncfiles_match(ind(dim_num(relevant).eq.dimsizes(variable_name)))
end if
end if
end if
return(ncfiles_wVar)
end;find_files_wVar
;===================================================================================================
; read in latlon variable from a file
; optionally interpolate onto a common 1x1 grid (to help with ensemble averaging)
undef("read_latlon_var")
function read_latlon_var(name_wfullpath:string,variable_name:string,REGRID:logical)
local name_wfullpath,variable_name,REGRID,fil
begin
;read in variable and perform some checks.
if .not.isfilepresent(name_wfullpath) then
;Hard exit if not correct, as only valid files should have been selected with find_files_wVar
print("There is not a appropriate file called "+name_wfullpath)
exit
end if
fil=addfile(name_wfullpath,"r")
if .not.isfilevar(fil,variable_name) then
;Hard exit if not correct, as only valid files should have been selected with find_files_wVar
print("There is not a "+variable_name+" variable in "+name_wfullpath)
exit
end if
var_in = fil->$variable_name$
;rename the (x,y) variables for later compatibility
if var_in!0.eq."latitude" then
var_in!0="lat"
end if
if var_in!1.eq."longitude" then
var_in!1="lon"
end if
if REGRID then
a = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")
lsdata = a->LSMASK
lat=lsdata&lat
lat@axis = "Y"
lat@long_name = "latitude"
lat@standard_name = "latitude"
lon=lsdata&lon
lon@axis = "X"
lon@long_name = "longitude"
lon@standard_name = "longitude"
var_out=new((/dimsizes(lat),dimsizes(lon)/),float,getFillValue(var_in))
var_out!0="lat"
var_out&lat=lat
var_out!1="lon"
var_out&lon=lon
copy_VarAtts(var_in,var_out)
var_out=linint2(var_in&lon,var_in&lat,var_in,True,lon,lat,0)
var_out@regrid="regridded from native grid using bilinear interpolation"
else
var_out=var_in
end if
return(var_out)
end;read_latlon_var
;===================================================================================================
; read in latlon variable from two files
; compute the difference in them
; optionally interpolate onto a common 1x1 grid (to help with ensemble averaging)
undef("read_diff_latlon_var")
function read_diff_latlon_var(name_wfullpath_a:string,name_wfullpath_b:string,variable_name:string,REGRID:logical)
; path for TS file(s), variable name, start year, and end year are read in.
local name_wfullpath_a,name_wfullpath_b,variable_name,REGRID,fil
begin
;read in variable and perform some checks.
if .not.isfilepresent(name_wfullpath_b).or..not.isfilepresent(name_wfullpath_a) then
;Hard exit if not correct, as only valid files should have been selected with find_files_wVar
print("There is either not a appropriate file called "+name_wfullpath_a+" or "+name_wfullpath_b)
exit
end if
fila=addfile(name_wfullpath_a,"r")
if .not.isfilevar(fila,variable_name) then
;Hard exit if not correct, as only valid files should have been selected with find_files_wVar
print("There is not a "+variable_name+" variable in "+name_wfullpath_a)
exit
end if
var_in_a = fila->$variable_name$
filb=addfile(name_wfullpath_b,"r")
if .not.isfilevar(filb,variable_name) then
;Hard exit if not correct, as only valid files should have been selected with find_files_wVar
print("There is not a "+variable_name+" variable in "+name_wfullpath_b)
exit
end if
var_in_b = filb->$variable_name$
diff_in=var_in_a
diff_in=var_in_a-var_in_b
diff_in@difference=name_wfullpath_a+" - "+name_wfullpath_b
;rename the (x,y) variables for later compatibility
if diff_in!0.eq."latitude" then
diff_in!0="lat"
end if
if diff_in!1.eq."longitude" then
diff_in!1="lon"
end if
if REGRID then
a = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")
lsdata = a->LSMASK
lat=lsdata&lat
lat@axis = "Y"
lat@long_name = "latitude"
lat@standard_name = "latitude"
lon=lsdata&lon
lon@axis = "X"
lon@long_name = "longitude"
lon@standard_name = "longitude"
diff_out=new((/dimsizes(lat),dimsizes(lon)/),float,getFillValue(diff_in))
diff_out!0="lat"
diff_out&lat=lat
diff_out!1="lon"
diff_out&lon=lon
copy_VarAtts(diff_in,diff_out)
diff_out=linint2(diff_in&lon,diff_in&lat,diff_in,True,lon,lat,0)
diff_out@regrid="regridded from native grid using bilinear interpolation"
else
diff_out=diff_in
end if
return(diff_out)
end;read_diff_latlon_var
;===================================================================================================
; extract a summary statistic for a region of a latlon variable from a file
undef("extract_latlon_areastat")
function extract_latlon_areastat(name_wfullpath:string,variable_name:string,lat_range:numeric,lon_range:numeric,stat_name:string)
local name_wfullpath,variable_name,fil,lat_range,lon_range,stat_name
begin
poss_stats=(/"mean","stddev"/)
if .not.any(poss_stats.eq.stat_name) then
print("extract_latlon_areastat: The only possible statistics currently accepted are "+str_join(poss_stats,","))
exit
end if
;read in variable and perform some checks.
if .not.isfilepresent(name_wfullpath) then
;Hard exit if not correct, as only valid files should have been selected with find_files_wVar
print("extract_latlon_areastat: There is not a appropriate file called "+name_wfullpath)
exit
end if
fil=addfile(name_wfullpath,"r")
if .not.isfilevar(fil,variable_name) then
;Hard exit if not correct, as only valid files should have been selected with find_files_wVar
print("extract_latlon_areastat: There is not a "+variable_name+" variable in "+name_wfullpath)
exit
end if
var_in = fil->$variable_name$
;rename the (x,y) variables for later compatibility
if var_in!0.eq."latitude" then
var_in!0="lat"
end if
if var_in!1.eq."longitude" then
var_in!1="lon"
end if
if .not.isdim(var_in,"lat") then
print("extract_latlon_areastat:"+variable_name+" does not have a latitude dimension in "+name_wfullpath)
exit
;not testing longitude as may need to apply this script for amoc at some point
end if
lon_range_posi=lon_range.ge.180.and.lon_range.le.360
lon_range_negi=lon_range.ge.-180.and.lon_range.lt.0
if any(lon_range_negi) then
if any(lon_range_posi) then
print("extract_latlon_areastat: Your lon_range is ("+lon_range(0)+","+lon_range(1)+") and I don't whether to use longitude of (-180:180) or (0:360)")
exit
else
var_in=lonFlip(var_in)
end if
end if
;extract subset and calculate lat. weights
wgts=NormCosWgtGlobe(var_in&lat)
wgts!0="lat"
wgts&lat=var_in&lat
subset=var_in({lat|lat_range(0):lat_range(1)},{lon|lon_range(0):lon_range(1)})
wgts_subset=wgts({lat_range(0):lat_range(1)})
;Calculate statistics
if stat_name.eq."mean".or.stat_name.eq."stddev" then
area_mn=wgt_areaave_Wrap(subset,wgts_subset,1.,0)
var_out=area_mn
if stat_name.eq."stddev" then
subset=(/(subset-area_mn)^2/)
var_out=sqrt(wgt_areaave(subset,wgts_subset,1.,0))
end if
end if
return(var_out)
end;read_latlon_var
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
undef("find_pair_files_wVar")
function find_pair_files_wVar(name_str_a:string,name_str_b:string,variable_name:string)
local name_str_a,name_str_b,variable_name,ncfile_wVar_a,ncfiles_wVar_b,match_count,name_i,ncfiles_wVar_both
begin
if isStrSubset(name_str_b,name_str_a).or.isStrSubset(name_str_a,name_str_b) then
print("find_pair_files_wVar: This program cannot work with one name string being a subset of the other")
print("find_pair_files_wVar: You may want to rename your cvdp_data files")
exit
end if
ncfiles_wVar_a=find_files_wVar(name_str_a,variable_name)
ncfiles_wVar_b=find_files_wVar(name_str_b,variable_name)
ncfiles_wVar_both=new((/dimsizes(ncfiles_wVar_a),2/),string) ;create a holding array for the matching ones
match_count=0
do name_i=0,dimsizes(ncfiles_wVar_a)-1
a_nopath=str_get_field(ncfiles_wVar_a(name_i),str_fields_count(ncfiles_wVar_a(name_i),"/"),"/")
gcm=str_get_field(str_sub_str(a_nopath,name_str_a,":"),1,":")
if .not.any(ismissing(str_match_ind(ncfiles_wVar_b,gcm))) then
if dimsizes(str_match_ind(ncfiles_wVar_b,gcm)).ne.1 then
print("find_pair_files_wVar: I've found multiple gcm id matches for "+gcm+" in "+name_str_a)
print("find_pair_files_wVar: I'm not sure how to continue.")
print("find_pair_files_wVar: matches = "+str_match(ncfiles_wVar_b,gcm))
exit
else
ncfiles_wVar_both(match_count,0)=ncfiles_wVar_a(name_i)
ncfiles_wVar_both(match_count,1)=str_match(ncfiles_wVar_b,gcm)
match_count=match_count+1
end if
end if
end do
if match_count.eq.0 then
print("find_pair_files_wVar: No matches found for cvdp_data names that differ solely by having "+name_str_b+" instead of "+name_str_a+" at the end")
ncfiles_wVar="missing"
else
ncfiles_wVar=ncfiles_wVar_both(ind(.not.ismissing(ncfiles_wVar_both(:,0))),:)
end if
return(ncfiles_wVar)
end;find_pair_files_wVar
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Program to run up a quick plot of a field
undef("roughlatlonplot")
procedure roughlatlonplot(Field)
local title,res,wks,plot
begin
if dimsizes(dimsizes(Field)).ne.2 then
print("Input only 2D latxlon arrays, please.")
return
end if
if isatt(Field,"long_name") then
title="rough plot of "+Field@long_name
else
title="rough plot"
end if
wks=gsn_open_wks("X11",title)
gsn_define_colormap(wks,"rainbow")
res = True ; plot mods desired
res@cnFillMode = "RasterFill" ; raster plot
res@cnFillOn = True
if isatt(Field,"cnLinesOn") then
res@cnLinesOn = Field@cnLinesOn
else
res@cnLinesOn = True
end if
res@cnLineLabelsOn = True
res@cnLevelSelectionMode = "AutomaticLevels"
res@gsnSpreadColors = True
res@gsnDraw = True
res@gsnFrame = True
res@gsnAddCyclic = False
if isatt(Field,"valid_range") then
res@cnLevelSelectionMode = "ManualLevels"
res@cnMaxLevelValF = Field@valid_range(1)
res@cnMinLevelValF = Field@valid_range(0)
res@cnLevelSpacingF = (Field@valid_range(1)-Field@valid_range(0))/10.
end if
if isatt(Field,"plotLatLons") then
res@mpLimitMode="LatLon"
res@mpMinLatF = Field@plotLatLons(0)
res@mpMaxLatF = Field@plotLatLons(1)
res@mpMinLonF = Field@plotLatLons(2)
res@mpMaxLonF = Field@plotLatLons(3)
end if
plot = gsn_csm_contour_map(wks,Field,res)
delete(wks)
end;roughlatlonplot
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;===================================================================================================
; read in a timeseries and return it along with its metadata (including @AnnCycle)
; options: ts_opt@renorm_climo=True alters when is the climatological period (set by ts_opt@renorm_climo_syear and ts_opt@renorm_climo_eyear)
; ts_opt@subset=True just strips out a subset of the timeseries [starting in ts_opt@subset_syear and ending ts_opt@subset_eyear)
; ts_opt@SEASON="DJF","MAM","JJA","SON" will convert a monthly timeseries into a seasonal one and return that. Also supports "ANN" for annual mean
; ts_opt@DETREND=True removes a linear trend from the timeseries
; ts_opt@RUN_STDDEV=True returns the running standard deviation of the timeseries rather than the timeseries itself
; It is caluclated over ts_opt@nstddev timepoints, which defaults to 30 years
; ts_opt@SMOOTH=True returns the running mean of the timeseries rather than the timeseries itself
; It is caluclated over ts_opt@nave timepoints, which defaults to 5 months for monthly fields, and 10 years for annual/seasonal ones
; ts_opt@make_absolute=True will add the Annual Cycle back onto the timeseries (stopping it being anomalies)
; ts_opt@NO_EXIT_ON_ERROR=True will override the default behaviour of failing upon error
undef("read_ts")
function read_ts(name_wfullpath:string,variable_name:string,ts_opt:logical)
local name_wfullpath,variable_name,ts_opt,fil,var_in,eyear,syear,file_years,file_syear,file_eyear,nmo,anncyc_correction,AnnCycle
begin
;read in variable and perform some checks.
if ts_opt.and.isatt(ts_opt,"NO_EXIT_ON_ERROR") then
NO_EXIT_ON_ERROR=ts_opt@NO_EXIT_ON_ERROR
missing=default_fillvalue("float")
missing@_FillValue=default_fillvalue("float")
else
NO_EXIT_ON_ERROR=False
end if
if .not.isfilepresent(name_wfullpath) then
;Hard exit if not correct, as only valid files should have been selected with find_files_wVar
if NO_EXIT_ON_ERROR then
return(missing)
else
print("There is not a appropriate file called "+name_wfullpath)
exit
end if
end if
fil=addfile(name_wfullpath,"r")
if any(variable_name.eq.(/"NWP","NCT"/)) then
if .not.isfilevar(fil,"nino34") then
if NO_EXIT_ON_ERROR then
return(missing)
else
;Hard exit if not correct, as only valid files should have been selected with find_files_wVar
print("There is not a "+variable_name+" variable in "+name_wfullpath)
exit
end if
end if
nino3=fil->nino3
nino4=fil->nino4
alpha=nino3*nino4 ;set up alpha array
alpha=where(alpha.gt.0,0.4,0)
var_in=nino3
if variable_name.eq."NWP" then
var_in=nino3-alpha*nino4;NCT
else
var_in=nino4-alpha*nino3;NWP
end if
delete([/nino3,nino4,alpha/])
else
if .not.isfilevar(fil,variable_name) then
if NO_EXIT_ON_ERROR then
return(missing)
else
;Hard exit if not correct, as only valid files should have been selected with find_files_wVar
print("There is not a "+variable_name+" variable in "+name_wfullpath)
exit
end if
end if
var_in = fil->$variable_name$
end if
read_ts_note=""
if .not.(isdim(var_in,"time").or.isdim(var_in,"TIME")) then
if NO_EXIT_ON_ERROR then
return(missing)
else
print("read_ts: "+variable_name+" does not have a time or TIME dimension. Suspect this is not a timeseries")
exit
end if
end if
if ts_opt.and.isatt(ts_opt,"make_absolute").and.ts_opt@make_absolute then
if .not.isStrSubset(variable_name,"monsoon") then
;ignore monsoon as they are already absolute variables
if .not.isatt(var_in,"AnnCycle").or.(dimsizes(var_in@AnnCycle).ne.12) then
if NO_EXIT_ON_ERROR then
return(missing)
else
print("read_ts: you have asked for make_absolute, but "+variable_name+" does not have a viable AnnCycle attribute")
exit
end if
else
AnnCycle=var_in@AnnCycle
do yr=0,dimsizes(var_in)-1,12
var_in(yr:yr+11) = (/ var_in(yr:yr+11)+AnnCycle/)
end do
end if
end if
end if
;Renormalised the climatology, if requested
if ts_opt.and.isatt(ts_opt,"renorm_climo").and.ts_opt@renorm_climo then
if .not.isatt(ts_opt,"renorm_climo_syear").or..not.isatt(ts_opt,"renorm_climo_eyear") then
print("read_ts: you have asked for renorm_climo, but not set the start or end year")
else
file_years=str_get_field(str_sub_str(name_wfullpath,"/../","/"),3,".")
file_syear = stringtointeger(str_get_field(file_years,1,"-"))
file_eyear = stringtointeger(str_get_field(file_years,2,"-"))
if (ts_opt@renorm_climo_syear.gt.0) then
check_custom_climo(name_wfullpath,file_syear,file_eyear,ts_opt@renorm_climo_syear,ts_opt@renorm_climo_eyear)
syear=ts_opt@renorm_climo_syear-file_syear
eyear=ts_opt@renorm_climo_eyear-file_syear
read_ts_note=read_ts_note+"Renormed climatology to "+ts_opt@renorm_climo_syear+"-"+ts_opt@renorm_climo_eyear+". "
else
;negative years are relative from end of the file
nyrs=(file_eyear-file_syear)+1
if abs(ts_opt@renorm_climo_syear).gt.nyrs then
print("read_ts: your relative climatolgy renormalisation starts before the file does")
else
syear=nyrs+ts_opt@renorm_climo_syear
eyear=nyrs+ts_opt@renorm_climo_eyear
read_ts_note=read_ts_note+"Renormed climatology to yrs "+(file_eyear+ts_opt@renorm_climo_syear)+"-"+(file_eyear+ts_opt@renorm_climo_eyear)+". "
end if
end if
if isdim(var_in,"time") then
;monthly timeseries
anncyc_correction = new(12,typeof(var_in),getFillValue(var_in))
do nmo=0,11
anncyc_correction(nmo) = dim_avg(var_in(nmo+(syear-1)*12:nmo+(eyear-1)*12:12))
end do
if isatt(var_in,"AnnCycle") then
AnnCycle=var_in@AnnCycle
AnnCycle=AnnCycle+anncyc_correction
var_in@AnnCycle=AnnCycle
end if
do yr=0,(file_eyear-file_syear),12
var_in(yr:yr+11) = (/ var_in(yr:yr+11)- anncyc_correction /)
end do
else
;annual timeseries
var_in=var_in-avg(var_in(syear-1:eyear))
end if
end if
end if
;take just a subset of the timeseries
if ts_opt.and.isatt(ts_opt,"subset").and.ts_opt@subset then
if .not.isatt(ts_opt,"subset_syear").or..not.isatt(ts_opt,"subset_eyear") then
print("read_ts: you have asked for 'subset', but not set start (subset_syear) or end (subset_eyear")
else
file_years=str_get_field(str_sub_str(name_wfullpath,"/../","/"),3,".")
file_syear = stringtointeger(str_get_field(file_years,1,"-"))
file_eyear = stringtointeger(str_get_field(file_years,2,"-"))
if (ts_opt@subset_syear.gt.0) then
check_custom_climo(name_wfullpath,file_syear,file_eyear,ts_opt@subset_syear,ts_opt@subset_eyear)
syear=ts_opt@subset_syear-file_syear
eyear=ts_opt@subset_eyear-file_syear
read_ts_note=read_ts_note+"Taken subset from "+ts_opt@subset_syear+" to "+ts_opt@subset_eyear+". "
else
;negative years are relative from end of the file
nyrs=(file_eyear-file_syear)+1
if abs(ts_opt@subset_syear).gt.nyrs then
print("read_ts: your relative temporal subset starts before the file does")
else
syear=nyrs+ts_opt@subset_syear
eyear=nyrs+ts_opt@subset_eyear
read_ts_note=read_ts_note+"Taken subset of yrs "+(file_eyear+ts_opt@subset_syear)+"-"+(file_eyear+ts_opt@subset_eyear)+". "
end if
end if
if isdim(var_in,"time") then ;monthly timeseries
var_out = var_in((syear-1)*12:eyear*12-1)
else ;annual timeseries
var_out=var_in(syear-1:eyear-1)
end if
end if
else
var_out=var_in ;return whole timeseries if subset not specified
end if
;Convert monthly to seasonal dataset
if ts_opt.and.isatt(ts_opt,"SEASON") then
if isdim(var_in,"time") then
;monthly timeseries
seasons=(/"MAM","JJA","SON","DJF","ANN"/)
if .not.any(seasons.eq.ts_opt@SEASON) then
print("read_ts: SEASON is set to convert to "+ts_opt@SEASON+" but I only take MAM, JJA, SON and DJF")
exit
else
if [email protected]."ANN" then
var_seas=var_out(5::12);set up array
var_seas=(/month_to_annual(var_out,1)/)
else
start_month=3*ind(seasons.eq.ts_opt@SEASON)-1
var_seas=var_out(start_month::12);set up array
var_seas=dim_avg_n((/var_out(start_month::12),var_out(start_month+1::12),var_out(start_month+2::12)/),0)
end if
delete(var_out)
var_out=var_seas;switch over variable
read_ts_note=read_ts_note+"Converted from monthly to "+ts_opt@SEASON+". "
end if
else
print("read_ts: SEASON is set to convert to "+ts_opt@SEASON+" but "+variable_name+" is not a monthly timseries")
exit
end if
end if
;remove a linear trend from the timeseries
if ts_opt.and.isatt(ts_opt,"DETREND").and.ts_opt@DETREND then
if isdim(var_out,"time") then
years=cd_calendar(var_out&time,4)
else
years=cd_calendar(var_out&TIME,4)
end if
var_out = dtrend_msg(years,var_out,True,False)
read_ts_note=read_ts_note+"A linear trend has been removed. "
end if
if ts_opt.and.isatt(ts_opt,"RUN_STDDEV").and.ts_opt@RUN_STDDEV then
if isatt(ts_opt,"nstddev") then
nstddev=ts_opt@nstddev
else
if isdim(var_out,"time").and.(.not.isatt(ts_opt,"SEASON")) then
nstddev=30*12 ;30 years
else
nstddev=30 ;30 years
end if
end if
ts_stddev=var_out
;set beginning and end to missing
ts_stddev(0:nstddev/2-1)=var_out@_FillValue
ts_stddev(dimsizes(ts_stddev)-nstddev/2:dimsizes(ts_stddev)-1)=var_out@_FillValue
do i=nstddev/2,dimsizes(ts_stddev)-nstddev/2-1
ts_stddev(i)=dim_stddev(var_out(i-nstddev/2:i+nstddev/2))
end do
var_out=(/ts_stddev/)
read_ts_note=read_ts_note+" The timeseries has been converted to a Running Standard Deviation over "+nstddev+" points. "
delete([/ts_stddev,nstddev/])
end if
if ts_opt.and.isatt(ts_opt,"SMOOTH").and.ts_opt@SMOOTH then
if isatt(ts_opt,"nave") then
nave=ts_opt@nave
else
if isdim(var_out,"time").and.(.not.isatt(ts_opt,"SEASON")) then
nave=5 ;5 month smoothing
else
nave=30 ;30 years
end if
end if
ts_smooth=var_out
;set beginning and end to missing
ts_smooth=runave(var_out, nave, 0)
var_out=(/ts_smooth/)
read_ts_note=read_ts_note+" The timeseries has been smoothed with a running mean over "+nave+" points. "
delete([/ts_smooth,nave/])
end if
var_out@read_ts_note=read_ts_note
return(var_out)
end;read_ts
;===================================================================================================
; return a single summary statistic for a timeseries
; stat_name = either one of "mean","stddev","skew","kurtosis" (as computed by stat4)
; or "pcvar", which gives the %-age variance explained (for EOF-based timeseries)
; or "AnnCycAmp" which gives max(AnnCycle)-min(AnnCycle)
undef("stat_ts_var")
function stat_ts_var(name_wfullpath:string,variable_name:string,stat_name_in:string,ts_opt:logical)
local name_wfullpath,variable_name,ts_opt,fil,poss_stats,stat_name
begin
;read in variable and perform some checks.
poss_stats=(/"mean","stddev","skew","kurtosis","pcvar","AnnCycAmp","stddev_bandpass_2-8yrs","dominant_period_interann","dominant_period"/)
if .not.any(poss_stats.eq.stat_name_in) then
print("stat_ts_var: The only possible statistics accepted are "+str_join(poss_stats,","))
exit
else
stat_name=stat_name_in
end if
var_in=read_ts(name_wfullpath,variable_name,ts_opt)
if stat_name.eq."stddev_bandpass_2-8yrs" then
; deploy a 2-8 yrs bandpass filter and set the variable to "stddev"
sigma = 1.0 ; Lanczos sigma
offset = 24 ; loose 2 years at each end
nwgt = 2*offset+1
low_thres = 1/(8.*12) ; freq for period of 8 yrs (1/months)
hi_thres = 1/(2.*12) ; freq for period of 2 yrs (1/months)
wgts_interann = filwgts_lanczos(nwgt,2,low_thres,hi_thres,sigma)
var_in=wgt_runave(var_in,wgts_interann,0)
stat_name="stddev"
end if
if .not.all(ismissing(var_in)) then
if (stat_name.eq."dominant_period_interann").or.(stat_name.eq."dominant_period") then
; determine number of datapoint per year
if isdim(var_in,"time") then
n_in_year=12.
else
n_in_year=1
end if
;calculate power spectra ;using most values from Brierley (2015, PlioMIP ENSO)
; but detrending will be done beforehand if requested, so iopt=0
sdof = specx_anal(var_in,0,5,0.10)
spectrum=(/sdof@spcx/n_in_year/)
frq=1./(sdof@frq*n_in_year)
if (stat_name.eq."dominant_period_interann") then
interann=ind(frq.ge.1.2.and.frq.lt.8)
interann_frq=frq(interann)
stat_out=interann_frq(maxind(spectrum(interann)))
else
stat_out=frq(maxind(spectrum))
end if
copy_VarAtts(var_in,stat_out)
stat_out@stat=stat_name
else
if stat_name.eq."AnnCycAmp" then
if .not.isatt(var_in,"AnnCycle") then
stat_out=default_fillvalue("float")
stat_out@_FillValue=default_fillvalue("float")
else
stat_out=max(var_in@AnnCycle)-min(var_in@AnnCycle)
end if
copy_VarAtts(var_in,stat_out)
stat_out@stat=stat_name
else
if stat_name.eq."pcvar" then
pcvar=str_squeeze(str_sub_str(var_in@pcvar,"%"," "))
stat_out=stringtofloat(pcvar)
else
stats=dim_stat4(var_in)
stat_ind=ind(poss_stats.eq.stat_name)
if stat_ind.gt.0.and.(all(var_in.eq.stats(0))) then
;If all the values are identical, only the mean is relevant
stat_out=default_fillvalue("float")
stat_out@_FillValue=default_fillvalue("float")
else
if stat_name.eq."stddev" then
stat_out=sqrt(stats(stat_ind))
else
stat_out=stats(stat_ind)
end if
end if
end if
copy_VarAtts(var_in,stat_out)
stat_out@stat=stat_name
end if
end if
else
stat_out=default_fillvalue("float")
stat_out@_FillValue=default_fillvalue("float")
end if
return(stat_out)
end;stat_ts_var
;===================================================================================================
; read in timeseries variable from multiple files
undef("read_ts_all")
function read_ts_all(name_str:string,variable_name:string,ts_opt:logical)
local name_str,variable_name,ncfiles,fil,var_in,file_years,file_syear,file_eyear
begin
;calculate the length of the longest timeseries (from the file and variable names)
ncfiles=find_files_wVar(name_str,variable_name)
max_years=0
names=ncfiles
do file_i=0,dimsizes(ncfiles)-1
file_years=str_get_field(str_sub_str(ncfiles(file_i),".cvdp_data.",":"),2,":")
file_syear = stringtointeger(str_get_field(file_years,1,"-."))
file_eyear = stringtointeger(str_get_field(file_years,2,"-."))
length_years=file_eyear-file_syear+1
if length_years.gt.max_years then
max_years=length_years
longest_file=file_i
end if
name_nopath=str_get_field(ncfiles(file_i),str_fields_count(ncfiles(file_i),"/"),"/")
names(file_i)=str_get_field(str_sub_str(name_nopath,".cvdp_data.",":"),1,":")
end do
;determine if field is monthly or annual
var_in = read_ts(ncfiles(longest_file),variable_name,ts_opt)
if isdim(var_in,"time") then
time=var_in&time
else
time=var_in&TIME
end if
var_out=new((/dimsizes(ncfiles),dimsizes(var_in)/),float)
var_out!0="runs"
var_out@run_names=names
var_out!1="time"
var_out&time=time
copy_VarAtts(var_in,var_out)
delete(var_in)
;read in the timeseries
do file_i=0,dimsizes(ncfiles)-1
var_in = read_ts(ncfiles(file_i),variable_name,ts_opt)
var_out(file_i,0:dimsizes(var_in)-1)=(/var_in/)
delete(var_in)
end do
return(var_out)
end;read_ts_all
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Plot obs vs ensemble mean maps ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
undef("plotCVDPcomparisonMaps")
procedure plotCVDPcomparisonMaps(strs,variable_name:string,wks_name,opt_res:logical,opt_panel_res:logical)
local mapres,panres,plots,ensmn,ensmn,patt,nvars,LeftStrings,RightStrings,CenterStrings,fnames,fnames,lat,lon,dims,opt_res,opt_panel_res
begin
; This has the supra-option of for opt_res NORMALISE and name_to_normalise_by
if opt_res.and.isatt(opt_res,"NORMALISE") then
;this option will plot the ratio of the fields
NORMALISE=opt_res@NORMALISE
delete(opt_res@NORMALISE)
if isatt(opt_res,"name_to_normalise_by") then
name_to_normalise_by=opt_res@name_to_normalise_by
delete(opt_res@name_to_normalise_by)
else
print("You have selected NORMALISE, but not provided the attribute name by which to normalise it (as name_to_normalise_by)")
exit
end if
else
NORMALISE=False
end if
;process some of the input variables
nvars=dimsizes(variable_name)
nstrs=dimsizes(strs)
plots=new(nstrs*nvars,graphic)
if isstring(wks_name) then
;Creating a new wks called wks_name
if isatt(wks_name,"filetype") then
wks=gsn_open_wks(wks_name@filetype,wks_name)
else
wks=gsn_open_wks("pdf",wks_name)
end if
else
;Adding this plot to pre-existing workstation
wks=wks_name
end if
mapres = True
;Altering some defaults
mapres@gsnDraw = False
mapres@gsnFrame = False
if nvars.eq.1 then
mapres@lbLabelBarOn = False
end if
mapres@cnLineLabelsOn = False
mapres@cnFillOn = True
mapres@cnLinesOn = False
mapres@cnFillPalette = "BlueDarkRed18"
mapres@mpGridAndLimbOn = True
mapres@mpGridLatSpacingF = 15. ; change latitude line spacing
mapres@mpGridLonSpacingF = 15. ; change longitude line spacing
;copy across the optional resources
copy_VarAtts(opt_res,mapres)
;using mpProjection as a major switch
if isatt(mapres,"mpProjection")[email protected]."WinkelTripel" then
mapres@vpYF = 0.95
mapres@vpHeightF = 0.3
mapres@vpXF = 0.2
mapres@vpWidthF = 0.6
mapres@gsnLeftStringOrthogonalPosF = -0.05
mapres@gsnLeftStringParallelPosF = .005
mapres@gsnRightStringOrthogonalPosF = -0.05
mapres@gsnRightStringParallelPosF = 0.96
mapres@gsnLeftStringFontHeightF = 0.02
mapres@gsnRightStringFontHeightF = 0.02
mapres@mpGeophysicalLineColor = "gray42"
mapres@mpPerimOn = False
mapres@mpGridLineColor = "transparent" ; trick ncl into drawing perimeter
end if
;load in the optional ones
if isatt(mapres,"gsnLeftString").and.dimsizes(mapres@gsnLeftString).eq.nstrs*nvars then
LeftStrings=mapres@gsnLeftString
delete(mapres@gsnLeftString)
end if
if isatt(mapres,"gsnRightString").and.dimsizes(mapres@gsnRightString).eq.nstrs*nvars then
RightStrings=mapres@gsnRightString
delete(mapres@gsnRightString)
end if
if isatt(mapres,"gsnCenterString").and.dimsizes(mapres@gsnCenterString).eq.nstrs*nvars then
CenterStrings=mapres@gsnCenterString
delete(mapres@gsnCenterString)
end if
if isatt(mapres,"tiMainString").and.dimsizes(mapres@tiMainString).eq.nstrs*nvars then
MainStrings=mapres@tiMainString
delete(mapres@tiMainString)
end if
panres=True
if opt_panel_res then
copy_VarAtts(opt_panel_res,panres)
end if
panres_options=(/"gsnMaximize","gsnPaperOrientation","gsnPanelYWhiteSpacePercent"/)
panres_new_defaults=True
panres_new_defaults@gsnMaximize = True
panres_new_defaults@gsnPaperOrientation = "portrait"
panres_new_defaults@gsnPanelYWhiteSpacePercent = 3.0
do i=0,dimsizes(panres_options)-1
if .not.isatt(panres,panres_options(i)) then
panres@$panres_options(i)$=panres_new_defaults@$panres_options(i)$
end if
end do
if nvars.eq.1 then
panres_singlevar_options=(/"gsnPanelLabelBar","pmLabelBarHeightF","pmLabelBarWidthF","lbBoxLineColor","lbTitleFontHeightF","lbTitlePosition","pmLabelBarOrthogonalPosF"/)
panres_new_defaults@gsnPanelLabelBar = True
panres_new_defaults@pmLabelBarHeightF = 0.05
panres_new_defaults@pmLabelBarWidthF = 0.55
panres_new_defaults@lbBoxLineColor = "gray70"
panres_new_defaults@lbTitleFontHeightF = 0.014
panres_new_defaults@lbTitlePosition = "Bottom"
panres_new_defaults@pmLabelBarOrthogonalPosF=-0.05
do i=0,dimsizes(panres_singlevar_options)-1
if .not.isatt(panres,panres_singlevar_options(i)) then
panres@$panres_singlevar_options(i)$=panres_new_defaults@$panres_singlevar_options(i)$
end if
end do
end if
;;determine names of files that have the correct field
do str_i=0,nstrs-1
do var_i=0,nvars-1
fnames=find_files_wVar(strs(str_i),variable_name(var_i))
dummy=read_latlon_var(fnames(0),variable_name(var_i),True)
lat=dummy&lat
lon=dummy&lon
ensmn=dummy
;load and average historical
patt=new((/dimsizes(fnames),dimsizes(lat),dimsizes(lon)/),float)
patt!0="run"
patt!1="lat"
patt&lat=lat
patt!2="lon"
patt&lon=lon
if NORMALISE then
norm_attrs=new(dimsizes(fnames),float)
end if
do mod_i=0,dimsizes(fnames)-1
if isatt(wks_name,"VERBOSE").and.(wks_name@VERBOSE) then
print("Loading "+variable_name(var_i)+" from file "+mod_i+":"+fnames(mod_i))
end if
if NORMALISE then
this_patt=read_latlon_var(fnames(mod_i),variable_name(var_i),True)
if (name_to_normalise_by.eq."relative_to_plot_area_average") then
norm_attrs(mod_i)=extract_latlon_areastat(fnames(mod_i),variable_name(var_i),(/mapres@mpMinLatF,mapres@mpMaxLatF/),(/mapres@mpMinLonF,mapres@mpMaxLonF/),"mean")
patt(mod_i,:,:)=(/this_patt-norm_attrs(mod_i)/)
else
norm_attrs(mod_i)=this_patt@$name_to_normalise_by$
patt(mod_i,:,:)=(/this_patt/norm_attrs(mod_i)/)
end if
delete(this_patt)
else
patt(mod_i,:,:)=read_latlon_var(fnames(mod_i),variable_name(var_i),True)
end if
end do
ensmn=(/dim_avg(patt(lat|:,lon|:,run|:))/)
;make the variable name the LabelBar title if not set
if .not.isatt(panres,"lbTitleString") then
panres@lbTitleString=dummy@long_name
if NORMALISE then
panres@lbTitleString="Normalised "+panres@lbTitleString
end if
end if
mapres@lbTitleString=panres@lbTitleString
;now create the ensemble mean plot
if isvar("LeftStrings") then
mapres@gsnLeftString = LeftStrings(nvars*str_i+var_i)
end if
if isvar("RightStrings") then
mapres@gsnRightString = RightStrings(nvars*str_i+var_i)
else
if NORMALISE then
mapres@gsnRightString = avg(norm_attrs)
end if
end if
if isvar("CenterStrings") then
mapres@gsnCenterString = CenterStrings(nvars*str_i+var_i)
end if
if isvar("MainStrings") then
mapres@tiMainString = MainStrings(nvars*str_i+var_i)
end if
plots(nvars*str_i+var_i)=gsn_csm_contour_map(wks,ensmn,mapres)
if NORMALISE then
delete([/patt,ensmn,dummy,fnames,norm_attrs/])
else
delete([/patt,ensmn,dummy,fnames/])
end if
end do;finish looping over the variables
end do
;now panel all the plots together...
if isatt(panres,"dims") then
dims=panres@dims
delete(panres@dims)
gsn_panel(wks,plots,dims,panres)
else
gsn_panel(wks,plots,(/nstrs,nvars/),panres)
end if
if isstring(wks_name) then