| 1 | ; |
|---|
| 2 | ; ---------------------- |
|---|
| 3 | ; Binning plot y=bin(x) |
|---|
| 4 | ; called by plt_map.pro |
|---|
| 5 | ; ---------------------- |
|---|
| 6 | ; |
|---|
| 7 | IF debug_w THEN print, " " |
|---|
| 8 | IF debug_w THEN print, " Enter ybinx..." |
|---|
| 9 | IF debug_w THEN print, " " |
|---|
| 10 | |
|---|
| 11 | ; organise bins: either read array or build form single number and min/max in fld_glo_mmx.def |
|---|
| 12 | ; fld2 is the field being bined |
|---|
| 13 | |
|---|
| 14 | IF (size(bin_interval))(0) EQ 1 THEN BEGIN |
|---|
| 15 | ; build... |
|---|
| 16 | ENDIF |
|---|
| 17 | |
|---|
| 18 | nbins = (size(bin_interval))(1) |
|---|
| 19 | |
|---|
| 20 | bin_plt = (bin_interval+shift(bin_interval, -1))/2 |
|---|
| 21 | bin_plt = [2*bin_plt(0)-bin_plt(1), bin_plt[0:nbins-2], 2*bin_plt(nbins-2)-bin_plt(nbins-3)] |
|---|
| 22 | |
|---|
| 23 | IF debug_w THEN print, ' bin_plt= ', bin_plt |
|---|
| 24 | IF debug_w THEN print, ' nbins= ', nbins |
|---|
| 25 | |
|---|
| 26 | ; 3rd field to bin regression ? |
|---|
| 27 | |
|---|
| 28 | IF var3_ybinx NE "" THEN sw3 = 1 ELSE sw3 = 0 |
|---|
| 29 | |
|---|
| 30 | ; mask fields with valmask for regression computation |
|---|
| 31 | |
|---|
| 32 | IF sw3 THEN BEGIN |
|---|
| 33 | idxmsk = where(fld GT valmask/10.) |
|---|
| 34 | IF idxmsk(0) NE -1 THEN fld(idxmsk) = valmask |
|---|
| 35 | IF idxmsk(0) NE -1 THEN fld3(idxmsk) = valmask |
|---|
| 36 | ENDIF |
|---|
| 37 | |
|---|
| 38 | ; ensure masks are the same |
|---|
| 39 | |
|---|
| 40 | idxmskpn = where(finite(fld, /nan)) |
|---|
| 41 | idxmskpn2 = where(finite(fld2, /nan)) |
|---|
| 42 | IF idxmskpn(0) NE -1 THEN fld(idxmskpn) = valmask |
|---|
| 43 | IF idxmskpn2(0) NE -1 THEN fld2(idxmskpn2) = valmask |
|---|
| 44 | |
|---|
| 45 | idxmskpn = where(fld GT valmask/10.) |
|---|
| 46 | idxmskpn2 = where(fld2 GT valmask/10.) |
|---|
| 47 | IF idxmskpn(0) NE -1 THEN fld2(idxmskpn) = valmask |
|---|
| 48 | IF idxmskpn2(0) NE -1 THEN fld(idxmskpn2) = valmask |
|---|
| 49 | |
|---|
| 50 | ; print min/max of field for debug |
|---|
| 51 | idxmskp = where(fld LE valmask/10.) |
|---|
| 52 | IF debug_w THEN print, ' Min/max fld= ', min(fld(idxmskp)), max(fld(idxmskp)) |
|---|
| 53 | IF debug_w THEN print, ' Min/max fld2= ', min(fld2(idxmskp)), max(fld2(idxmskp)) |
|---|
| 54 | |
|---|
| 55 | ; remove mean seasonal cycle if required |
|---|
| 56 | |
|---|
| 57 | IF cmd.trend GT 0 THEN BEGIN |
|---|
| 58 | fld = trends(fld, cmd.trend, 'xyt') |
|---|
| 59 | IF sw3 THEN fld3 = trends(fld3, cmd.trend, 'xyt') |
|---|
| 60 | ENDIF |
|---|
| 61 | IF cmd2.trend GT 0 THEN BEGIN |
|---|
| 62 | fld2 = trends(fld2, cmd2.trend, 'xyt') |
|---|
| 63 | ENDIF |
|---|
| 64 | |
|---|
| 65 | ; select months if required |
|---|
| 66 | |
|---|
| 67 | CASE cmd.timave OF |
|---|
| 68 | '1d': ntxt = "All days" |
|---|
| 69 | ELSE :ntxt = "All months" |
|---|
| 70 | ENDCASE |
|---|
| 71 | |
|---|
| 72 | IF stddev_mth NE '00' THEN BEGIN |
|---|
| 73 | @mth_decode |
|---|
| 74 | |
|---|
| 75 | fldn = (fld)[*, *, reform(idxm(0,*), njpt)] |
|---|
| 76 | fld2n = (fld2)[*, *, reform(idxm(0,*), njpt)] |
|---|
| 77 | IF sw3 THEN fld3n = (fld3)[*, *, reform(idxm(0,*), njpt)] |
|---|
| 78 | |
|---|
| 79 | FOR imth = 1, nmth-1 DO BEGIN |
|---|
| 80 | |
|---|
| 81 | fldn = [fldn, (fld)[*, *, reform(idxm(imth,*), njpt)]] |
|---|
| 82 | fld2n = [fld2n, (fld2)[*, *, reform(idxm(imth,*), njpt)]] |
|---|
| 83 | IF sw3 THEN fld3n = [fld3n, (fld3)[*, *, reform(idxm(imth,*), njpt)]] |
|---|
| 84 | |
|---|
| 85 | ENDFOR |
|---|
| 86 | jpt = njpt |
|---|
| 87 | |
|---|
| 88 | fld = fldn |
|---|
| 89 | fld2 = fld2n |
|---|
| 90 | IF sw3 THEN fld3 = fld3n |
|---|
| 91 | ENDIF |
|---|
| 92 | |
|---|
| 93 | ; for now just 2d fields |
|---|
| 94 | ; IF nzt NE 1 THEN BEGIN |
|---|
| 95 | ; print, '***** 2D field only for now in ybinx ****' |
|---|
| 96 | ; stop |
|---|
| 97 | ; ENDIF |
|---|
| 98 | |
|---|
| 99 | ; find indexes of var2 in each bin |
|---|
| 100 | |
|---|
| 101 | idxb = lonarr(nbins+1, nxt*nyt*jpt) |
|---|
| 102 | binpop = lonarr(nbins+1) |
|---|
| 103 | |
|---|
| 104 | ib = 1 |
|---|
| 105 | |
|---|
| 106 | WHILE ib LE nbins-1 DO BEGIN |
|---|
| 107 | |
|---|
| 108 | indices = where(fld2 GT bin_interval(ib-1) AND fld2 LE bin_interval(ib)) |
|---|
| 109 | binpop(ib) = n_elements(indices) |
|---|
| 110 | idxb[ib, 0:binpop(ib)-1] = where(fld2 GT bin_interval(ib-1) AND fld2 LE bin_interval(ib)) |
|---|
| 111 | IF debug_w THEN print, ' Size of bin(i) ', ib, binpop(ib) |
|---|
| 112 | |
|---|
| 113 | ib = ib + 1 |
|---|
| 114 | ENDWHILE |
|---|
| 115 | |
|---|
| 116 | index0 = where(fld2 LE bin_interval(0)) |
|---|
| 117 | binpop(0) = n_elements(index0) |
|---|
| 118 | idxb(0, 0:binpop(0)-1) = index0 |
|---|
| 119 | |
|---|
| 120 | indexm = n_elements(where(fld2 GT bin_interval(nbins-1))) |
|---|
| 121 | binpop(nbins) = n_elements(indexm) |
|---|
| 122 | idxb(nbins, 0:binpop(nbins)-1) = indexm |
|---|
| 123 | |
|---|
| 124 | ; compute maximum number of indices for one bin |
|---|
| 125 | |
|---|
| 126 | max_idx = max(binpop) |
|---|
| 127 | |
|---|
| 128 | IF debug_w THEN print, ' Max bin size ', max_idx |
|---|
| 129 | IF debug_w THEN print, ' Number of bins ', nbins |
|---|
| 130 | fldy = fltarr(nbins+1, max_idx) |
|---|
| 131 | fldy[*, *] = !values.f_nan |
|---|
| 132 | fldys = fltarr(nbins+1, max_idx) |
|---|
| 133 | fldys[*, *] = !values.f_nan |
|---|
| 134 | surfb = fltarr(nbins+1, max_idx) |
|---|
| 135 | surfb[*, *] = !values.f_nan |
|---|
| 136 | |
|---|
| 137 | IF sw3 THEN BEGIN |
|---|
| 138 | fldy2 = fltarr(nbins+1, max_idx) |
|---|
| 139 | fldy2[*, *] = !values.f_nan |
|---|
| 140 | ENDIF |
|---|
| 141 | |
|---|
| 142 | ; extract fld1/fld3 arrays in each bin |
|---|
| 143 | |
|---|
| 144 | ; first build fld*e1te2t for later average |
|---|
| 145 | |
|---|
| 146 | surf = reform(e1t(firstxt:lastxt,firstyt:lastyt)*e2t(firstxt:lastxt,firstyt:lastyt), nxt*nyt) |
|---|
| 147 | surf = reform(surf#replicate(1, jpt), nxt, nyt, jpt) |
|---|
| 148 | |
|---|
| 149 | flds = fld*surf |
|---|
| 150 | |
|---|
| 151 | ib = 0 |
|---|
| 152 | WHILE ib LE nbins DO BEGIN |
|---|
| 153 | |
|---|
| 154 | IF debug_w THEN print, 'bin = ', ib |
|---|
| 155 | binsz = binpop(ib) |
|---|
| 156 | IF binsz GT 1 THEN BEGIN |
|---|
| 157 | fldy(ib, 0:binsz-1) = fld(idxb(ib, 0:binsz-1)) |
|---|
| 158 | IF debug_w THEN print, 'fld(idxb(ib, 0:binsz-1)) =',fld(idxb(ib, 0:binsz-1)) |
|---|
| 159 | fldys(ib, 0:binsz-1) = flds(idxb(ib, 0:binsz-1)) |
|---|
| 160 | surfb(ib, 0:binsz-1) = surf(idxb(ib, 0:binsz-1)) |
|---|
| 161 | IF sw3 THEN fldy2(ib, 0:binsz-1) = fld3(idxb(ib, 0:binsz-1)) |
|---|
| 162 | ENDIF |
|---|
| 163 | ib = ib + 1 |
|---|
| 164 | ENDWHILE |
|---|
| 165 | |
|---|
| 166 | |
|---|
| 167 | yplt = fltarr(nbins+1) |
|---|
| 168 | yerr = fltarr(nbins+1) |
|---|
| 169 | |
|---|
| 170 | IF NOT sw3 THEN BEGIN |
|---|
| 171 | |
|---|
| 172 | ; if bining only fld1, compute average in each bin and plot |
|---|
| 173 | ; average |
|---|
| 174 | |
|---|
| 175 | mean_fld = 0. |
|---|
| 176 | |
|---|
| 177 | ib = 0 |
|---|
| 178 | WHILE ib LE nbins DO BEGIN |
|---|
| 179 | |
|---|
| 180 | binsz = binpop(ib) |
|---|
| 181 | |
|---|
| 182 | IF binsz GT 1 THEN BEGIN |
|---|
| 183 | |
|---|
| 184 | sfc_tot = total(surfb(ib, 0:binsz-1)) |
|---|
| 185 | yplt(ib) = total(fldys(ib, 0:binsz-1))/sfc_tot |
|---|
| 186 | yerr(ib) = sqrt((moment(fldy(ib, 0:binsz-1)))[1]) |
|---|
| 187 | |
|---|
| 188 | mean_fld = mean_fld + yplt(ib)*float(binpop(ib)) |
|---|
| 189 | |
|---|
| 190 | ; print bin info |
|---|
| 191 | IF ib GT 0 AND ib LT nbins THEN print, ' Bin size, occurence, average: ', bin_interval(ib-1), bin_interval(ib), binpop(ib), (binpop(ib)/total(binpop))*100., yplt(ib) |
|---|
| 192 | IF ib EQ 0 THEN print, ' Bin size, occurence, average: min' , bin_interval(ib), binpop(ib), (binpop(ib)/total(binpop))*100. , yplt(ib) |
|---|
| 193 | IF ib EQ nbins THEN print, ' Bin size, occurence, average: ', bin_interval(ib-1),' max ', binpop(ib), (binpop(ib)/total(binpop))*100. , yplt(ib) |
|---|
| 194 | ENDIF ELSE yplt(ib) = !values.f_nan |
|---|
| 195 | |
|---|
| 196 | ib = ib + 1 |
|---|
| 197 | ENDWHILE |
|---|
| 198 | |
|---|
| 199 | mean_fld = mean_fld/total(binpop) |
|---|
| 200 | |
|---|
| 201 | print, ' Total number of bins = ', total(binpop) |
|---|
| 202 | print, ' Mean field in domain = ', mean_fld |
|---|
| 203 | |
|---|
| 204 | ; specify plot attributes |
|---|
| 205 | varname = varlegend2 |
|---|
| 206 | |
|---|
| 207 | ENDIF ELSE BEGIN |
|---|
| 208 | ; if bining fld1=f(fld3), compute regression in each bin and plot |
|---|
| 209 | |
|---|
| 210 | ; compute regression for each bin |
|---|
| 211 | |
|---|
| 212 | mean_fld = 0. |
|---|
| 213 | |
|---|
| 214 | ib = 0 |
|---|
| 215 | WHILE ib LE nbins DO BEGIN |
|---|
| 216 | |
|---|
| 217 | binsz = binpop(ib) |
|---|
| 218 | |
|---|
| 219 | IF binsz GT 1 THEN BEGIN |
|---|
| 220 | |
|---|
| 221 | idx1 = where(fldy(ib, 0:binsz-1) NE valmask) |
|---|
| 222 | idx2 = where(fldy2(ib, 0:binsz-1) NE valmask) |
|---|
| 223 | tab1 = (fldy(ib, 0:binsz-1))(idx1) |
|---|
| 224 | tab2 = (fldy2(ib, 0:binsz-1))(idx2) |
|---|
| 225 | coeff = linfit(tab2, tab1, CHISQ = linerr, PROB = proberr, SIGMA = sigmaerr) |
|---|
| 226 | yplt(ib) = coeff(1) |
|---|
| 227 | yerr(ib) = sigmaerr(1) |
|---|
| 228 | |
|---|
| 229 | mean_fld = mean_fld + yplt(ib)*float(binpop(ib)) |
|---|
| 230 | |
|---|
| 231 | ; print bin info |
|---|
| 232 | IF ib GT 0 AND ib LT nbins THEN print, ' Bin size, occurence, regress.: ', bin_interval(ib-1), bin_interval(ib), binpop(ib), (binpop(ib)/total(binpop))*100., yplt(ib) |
|---|
| 233 | IF ib EQ 0 THEN print, ' Bin size, occurence, regress.: min' , bin_interval(ib), binpop(ib), (binpop(ib)/total(binpop))*100. , yplt(ib) |
|---|
| 234 | IF ib EQ nbins THEN print, ' Bin size, occurence, regress.: ', bin_interval(ib-1),' max ', binpop(ib), (binpop(ib)/total(binpop))*100. , yplt(ib) |
|---|
| 235 | |
|---|
| 236 | ENDIF ELSE yplt(ib) = !values.f_nan |
|---|
| 237 | |
|---|
| 238 | ib = ib + 1 |
|---|
| 239 | ENDWHILE |
|---|
| 240 | |
|---|
| 241 | mean_fld = mean_fld/total(binpop) |
|---|
| 242 | |
|---|
| 243 | print, ' Total number of bins = ', total(binpop) |
|---|
| 244 | print, ' Mean regression in domain = ', mean_fld |
|---|
| 245 | |
|---|
| 246 | ; specify plot attributes |
|---|
| 247 | varname = 'Regression of '+field.name+' and '+field.name3 |
|---|
| 248 | |
|---|
| 249 | ENDELSE |
|---|
| 250 | |
|---|
| 251 | ; plot |
|---|
| 252 | |
|---|
| 253 | ; define text, line color, thickness and type |
|---|
| 254 | |
|---|
| 255 | boxybinx = def_box(cmd.plt, dimplot, leg_name, time_stride) |
|---|
| 256 | vardate = date_txt |
|---|
| 257 | varunit = ' ['+ntxt+' in box '+leg_name+']' |
|---|
| 258 | |
|---|
| 259 | overc = overlay_type(iover, dimplot) |
|---|
| 260 | |
|---|
| 261 | pltcmd = 'pltsc,yplt,bin_plt,minc,maxc,minc2,maxc2,varlegend'+com_strplt+overc+',STY1D=-6,subtitle=""' |
|---|
| 262 | |
|---|
| 263 | printf, nulhis, strcompress(pltcmd) |
|---|
| 264 | IF debug_w THEN print, ' ', pltcmd |
|---|
| 265 | res = execute(pltcmd) |
|---|
| 266 | |
|---|
| 267 | ; plot +/- 1 stdedv for field binnig |
|---|
| 268 | |
|---|
| 269 | pltcmd = 'pltsc,yplt-yerr,bin_plt,minc,maxc,minc2,maxc2,varlegend'+com_strplt+',ov1d=1,COLOR=1, thick=1, STY1D=-1,subtitle=""' |
|---|
| 270 | res = execute(pltcmd) |
|---|
| 271 | pltcmd = 'pltsc,yplt+yerr,bin_plt,minc,maxc,minc2,maxc2,varlegend'+com_strplt+',ov1d=1,COLOR=1, thick=1, STY1D=-1,subtitle=""' |
|---|
| 272 | res = execute(pltcmd) |
|---|
| 273 | |
|---|
| 274 | IF debug_w THEN print, " ... Exit ybinx" |
|---|
