-
Notifications
You must be signed in to change notification settings - Fork 0
/
a_ond_g5ec_ens_mems_fc_cat_probs.R
137 lines (101 loc) · 3.94 KB
/
a_ond_g5ec_ens_mems_fc_cat_probs.R
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
# R code to read NMME ensemble members stored in a netcdf file downloaded from
# IRI data library and calculate ensemble probabilities
# lats, lons, years set by expert script
#
# probabilities stored in cpt format for imput into cpt probability forecast verification (PFV)
#
# install.packages(c("ncdf4"))
library(ncdf4)
tercile = TRUE # quintile if FALSE
glosea5 = TRUE # ECMWF if FALSE
if (glosea5 & tercile) probs_file='g5_ppn_ke_aug_OND2020_etprobs.txt'
if (glosea5 & (tercile == FALSE)) probs_file='g5_ppn_ke_aug_OND2020_eqprobs.txt'
if ((glosea5 == FALSE) & tercile) probs_file='ec_ppn_ke_aug_OND2020_etprobs.txt'
if ((glosea5 == FALSE) & (tercile == FALSE)) probs_file='ec_ppn_ke_aug_OND2020_eqprobs.txt'
tmths= "-01-01"
# nem=10
if (glosea5) hc_name <- "g5_ppn_ke_aug_OND9316_mems"
if (glosea5 == FALSE) hc_name <- "ec_ppn_ke_aug_OND9316_mems"
hc_fname <- paste(hc_name, ".nc", sep = "")
# open NetCDF file containing NMME ensemble hindcast members
nchc <- nc_open(hc_fname)
print(nchc)
lons <- ncvar_get(nchc, "X")
nlon <- dim(lons)
lats <- ncvar_get(nchc, "Y")
nlat <- dim(lats)
yrs <- ncvar_get(nchc, "S")
ny <- dim(yrs)
mems <- ncvar_get(nchc,"M")
nem <- dim(mems)
if (glosea5) prec0 <- ncvar_get(nchc, "precipitation_flux")
if (glosea5 == FALSE) prec0 <- ncvar_get(nchc, "prcp")
p1 <- array(prec0, , dim=c(nlon,nlat,nem,ny))
if (glosea5) fc_name <- "g5_ppn_ke_aug_OND2020_mems"
if (glosea5 == FALSE) fc_name <- "ec_ppn_ke_aug_OND2020_mems"
fc_fname <- paste(fc_name, ".nc", sep = "")
# open NetCDF file containing NMME ensemble hindcast members
ncfc <- nc_open(fc_fname)
print(ncfc)
fmems <- ncvar_get(ncfc,"M")
nemf <- dim(fmems)
fyrs <- ncvar_get(ncfc, "S")
nyf <- dim(fyrs)
if (glosea5) precf <- ncvar_get(ncfc, "precipitation_flux")
if (glosea5 == FALSE) precf <- ncvar_get(ncfc, "prcp")
pf <- array(precf, , dim=c(nlon,nlat,nemf))
# Loop to evaluate requency probs
probs <- array(0,, dim=c(nlon,nlat,3))
for (lon1 in 1:nlon) {
for (lat1 in 1:nlat) {
p2=p1[lon1,lat1,,]
p2f=pf[lon1,lat1,]
if (tercile) tce = quantile(p2, c(0.333,0.667))
if (tercile == FALSE) tce = quantile(p2, c(0.2,0.8))
p3= (p2f > tce[1]) + (p2f > tce[2])
probs[lon1,lat1,1] = 100*sum(p3 == 0)/nemf
probs[lon1,lat1,3] = 100*sum(p3 == 2)/nemf
probs[lon1,lat1,2]= 100.0-(probs[lon1,lat1,1]+probs[lon1,lat1,3])
}
}
# SAVE PROBS in CPT format
lines = c()
lines = c(lines, paste("xmlns:cpt=http://iri.columbia.edu/CPT/v10/"))
lines = c(lines, paste("cpt:ncats=3"))
if (tercile) h1 = "cpt:field=precipitation, cpt:C=1, cpt:clim_prob=0.333333333333, cpt:T="
if (tercile == FALSE) h1 = "cpt:field=precipitation, cpt:C=1, cpt:clim_prob=0.2, cpt:T="
h1 = paste0(h1, toString(as.integer(fyrs/12+1960)))
h1 = paste0(h1, tmths)
h1 = paste(h1, ", cpt:nrow=")
h1 = paste0(h1, toString(nlat))
h1 = paste(h1, ", cpt:ncol=")
h1 = paste0(h1, toString(nlon))
h1 = paste(h1, ", cpt:row=Y, cpt:col=X, cpt:units=probability (%), cpt:missing=-1.00000000000")
h2 = paste(lons,collapse=" ")
lines = c(lines, h1)
lines =c(lines, h2)
for (lat1 in 1:nlat) {
h3= paste(probs[,lat1,1],collapse=" ")
h3= paste(lats[lat1], h3)
lines = c(lines,h3)
}
if (tercile) lines=c(lines, "cpt:C=2, cpt:clim_prob=0.333333333334")
if (tercile == FALSE) lines=c(lines, "cpt:C=2, cpt:clim_prob=0.6")
lines =c(lines, h2)
for (lat1 in 1:nlat) {
h3= paste(probs[,lat1,2],collapse=" ")
h3= paste(lats[lat1], h3)
lines = c(lines,h3)
}
if (tercile) lines=c(lines, "cpt:C=3, cpt:clim_prob=0.333333333333")
if (tercile == FALSE) lines=c(lines, "cpt:C=3, cpt:clim_prob=0.2")
lines =c(lines, h2)
for (lat1 in 1:nlat) {
h3= paste(probs[,lat1,3],collapse=" ")
h3= paste(lats[lat1], h3)
lines = c(lines,h3)
}
# open file, write lines, close file
out = file(probs_file)
writeLines(lines, out)
close(out)