clamp: lai/99.ncl@0c6405ab2ff4

     1 ;********************************************************

     2 ; histogram normalized by rain and compute correleration

     3 ;********************************************************

     4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

     5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

     6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

     8 procedure pminmax(data:numeric,name:string)

     9 begin

    10   print ("min/max " + name + " = " + min(data) + "/" + max(data))

    11   if(isatt(data,"units")) then

    12     print (name + " units = " + data@units)

    13   end if

    14 end

    16 ; Main code.

    17 begin

    19 ;nclass = 18

    20  nclass = 20

    22  plot_type     = "ps"

    23  plot_type_new = "png"

    25 ;************************************************

    26 ; read in data: model

    27 ;************************************************

    29  model_name = "b30.061n"

    30  model_grid = "T31"

    32  dirm  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    33  film  = "b30.061n_1995-2004_MONS_climo_lnd.nc"

    34 ;film  = "i01.03cn_1545-1569_MONS_climo.nc"

    35  fm    = addfile(dirm+film,"r")

    37  laimod  = fm->TLAI

    39 ;************************************************

    40 ; read in data: observed

    41 ;************************************************

    43  ob_name = "MODIS MOD 15A2 2000-2005"

    45  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"

    46  filo1  = "land_class_"+model_grid+".nc"

    47  filo2  = "LAI_2000-2005_MONS_"+model_grid+".nc"

    49  fo1 = addfile(diro+filo1,"r")

    50  fo2 = addfile(diro+filo2,"r")

    52  classob    = tofloat(fo1->LAND_CLASS)

    53  laiob      = fo2->LAI

    54 ;*******************************************************************

    55 ; Calculate "nice" bins for binning the data in equally spaced ranges

    56 ;********************************************************************

    57   nclassn     = nclass + 1

    58   range       = fspan(0,nclassn-1,nclassn)

    59 ; print (range)

    61 ; Use this range information to grab all the values in a

    62 ; particular range, and then take an average.

    64   nr           = dimsizes(range)

    65   nx           = nr-1

    66   xvalues      = new((/2,nx/),float)

    67   xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.

    68   dx           = xvalues(0,1) - xvalues(0,0)       ; range width

    69   dx4          = dx/4                              ; 1/4 of the range

    70   xvalues(1,:) = xvalues(0,:) - dx/5.

    71 ;-----------------------------------------------------------------

    73 ;-----------------------------------------------------------------

    74 ;(B) max

    75 ;--------------------------------------------------------------------

    76 ; get data

    78 ; observed

    79   laiob_max = laiob(0,:,:)

    80   s         = laiob(:,0,0)

    81   laiob_max@long_name = "Leaf Area Index Max"

    83   dsizes_z = dimsizes(laiob)

    84   nlat     = dsizes_z(1)

    85   nlon     = dsizes_z(2)

    87   do j = 0,nlat-1

    88   do i = 0,nlon-1

    89      s = laiob(:,j,i)

    90      laiob_max(j,i) = max(s)

    91   end do

    92   end do

    94 ; print (min(y)+"/"+max(y))

    95   delete (s)

    96   delete (dsizes_z)

    97 ;-------------------------

    98 ; model

    99   laimod_max = laimod(0,:,:)

   100   s          = laimod(:,0,0)

   101   laimod_max@long_name = "Leaf Area Index Max"

   103   dsizes_z = dimsizes(laimod)

   104   nlat     = dsizes_z(1)

   105   nlon     = dsizes_z(2)

   107   do j = 0,nlat-1

   108   do i = 0,nlon-1

   109      s = laimod(:,j,i)

   110      laimod_max(j,i) = max(s)

   111   end do

   112   end do

   114 ; print (min(laimod_max)+"/"+max(laimod_max))

   115   delete (s)

   116   delete (dsizes_z)

   117 ;------------------------

   118   DATA11_1D = ndtooned(classob)

   119   DATA12_1D = ndtooned(laiob_max)

   120   DATA22_1D = ndtooned(laimod_max)

   122   yvalues      = new((/2,nx/),float)

   123   mn_yvalues   = new((/2,nx/),float)

   124   mx_yvalues   = new((/2,nx/),float)

   126   do nd=0,1

   128 ; See if we are doing model or observational data.

   130     if(nd.eq.0) then

   131       data_ob  = DATA11_1D

   132       data_mod = DATA12_1D

   133     else

   134       data_ob  = DATA11_1D

   135       data_mod = DATA22_1D

   136     end if

   138 ; Loop through each range and check for values.

   140     do i=0,nr-2

   141       if (i.ne.(nr-2)) then

   142 ;        print("")

   143 ;        print("In range ["+range(i)+","+range(i+1)+")")

   144         idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))

   145       else

   146 ;        print("")

   147 ;        print("In range ["+range(i)+",)")

   148         idx = ind(range(i).le.data_ob)

   149       end if

   151 ; Calculate average, and get min and max.

   153       if(.not.any(ismissing(idx))) then

   154         yvalues(nd,i)    = avg(data_mod(idx))

   155         mn_yvalues(nd,i) = min(data_mod(idx))

   156         mx_yvalues(nd,i) = max(data_mod(idx))

   157         count = dimsizes(idx)

   158       else

   159         count            = 0

   160         yvalues(nd,i)    = yvalues@_FillValue

   161         mn_yvalues(nd,i) = yvalues@_FillValue

   162         mx_yvalues(nd,i) = yvalues@_FillValue

   163       end if

   165 ;     print(nd + ": " + count + " points, avg = " + yvalues(nd,i))

   166 ;     print("Min/Max:  " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))

   168 ; Clean up for next time in loop.

   170       delete(idx)

   171     end do

   172     delete(data_ob)

   173     delete(data_mod)

   174   end do

   175 ;-----------------------------------------------------------------

   176 ; compute correlation coef and M score

   178   u = yvalues(0,:)

   179   v = yvalues(1,:)

   181   good = ind(.not.ismissing(u) .and. .not.ismissing(v))

   182   uu = u(good)

   183   vv = v(good)

   185   ccrMax = esccr(uu,vv,0)

   186 ; print (ccrMax)

   188 ; new eq

   189   bias = sum(abs(vv-uu)/(vv+uu))

   190   Mmax = (1.- (bias/dimsizes(uu)))*5.

   192   print (Mmax)

   194  delete (u)

   195  delete (v)

   196  delete (uu)

   197  delete (vv)

   198 ;--------------------------------------------------------------------

   199 ; histogram res

   201   resm                = True

   202   resm@gsnMaximize    = True

   203   resm@gsnDraw        = False

   204   resm@gsnFrame       = False

   205   resm@xyMarkLineMode = "Markers"

   206   resm@xyMarkerSizeF  = 0.014

   207   resm@xyMarker       = 16

   208   resm@xyMarkerColors = (/"Brown","Blue"/)

   209 ; resm@trYMinF        = min(mn_yvalues) - 10.

   210 ; resm@trYMaxF        = max(mx_yvalues) + 10.

   211   resm@trYMinF        = min(mn_yvalues) - 2

   212   resm@trYMaxF        = max(mx_yvalues) + 4

   214   resm@tiYAxisString  = "Max LAI (Leaf Area Index)"

   215   resm@tiXAxisString  = "Land Cover Type"

   217   max_bar = new((/2,nx/),graphic)

   218   min_bar = new((/2,nx/),graphic)

   219   max_cap = new((/2,nx/),graphic)

   220   min_cap = new((/2,nx/),graphic)

   222   lnres = True

   223   line_colors = (/"brown","blue"/)

   224 ;------------------------------------------------------------------

   225 ; Start the graphics.

   227   plot_name = "histogram_max"

   228   title     = model_name + " vs Observed"

   229   resm@tiMainString  = title

   231   wks   = gsn_open_wks (plot_type,plot_name)

   232 ;-----------------------------

   233 ; Add a boxed legend using the more simple method

   235   resm@pmLegendDisplayMode    = "Always"

   236 ; resm@pmLegendWidthF         = 0.1

   237   resm@pmLegendWidthF         = 0.08

   238   resm@pmLegendHeightF        = 0.05

   239   resm@pmLegendOrthogonalPosF = -1.17

   240 ; resm@pmLegendOrthogonalPosF = -1.00  ;(downward)

   241 ; resm@pmLegendParallelPosF   =  0.18

   242   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

   244 ; resm@lgPerimOn              = False

   245   resm@lgLabelFontHeightF     = 0.015

   246   resm@xyExplicitLegendLabels = (/"observed",model_name/)

   247 ;-----------------------------

   248   tRes  = True

   249   tRes@txFontHeightF = 0.025

   251   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"

   253   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

   255   xy = gsn_csm_xy(wks,xvalues,yvalues,resm)

   256 ;-------------------------------

   257 ;Attach the vertical bar and the horizontal cap line

   259   do nd=0,1

   260     lnres@gsLineColor = line_colors(nd)

   261     do i=0,nx-1

   263       if(.not.ismissing(mn_yvalues(nd,i)).and. \

   264          .not.ismissing(mx_yvalues(nd,i))) then

   266 ; Attach the vertical bar, both above and below the marker.

   268         x1 = xvalues(nd,i)

   269         y1 = yvalues(nd,i)

   270         y2 = mn_yvalues(nd,i)

   271         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

   273         y2 = mx_yvalues(nd,i)

   274         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

   276 ; Attach the horizontal cap line, both above and below the marker.

   278         x1 = xvalues(nd,i) - dx4

   279         x2 = xvalues(nd,i) + dx4

   280         y1 = mn_yvalues(nd,i)

   281         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

   283         y1 = mx_yvalues(nd,i)

   284         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

   285       end if

   286     end do

   287   end do

   289   draw(xy)

   290   frame(wks)

   292   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   293 ; system("rm "+plot_name+"."+plot_type)

   294 ; system("rm "+plot_name+"-1."+plot_type_new)

   295 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   297   clear (wks)

   299  delete (DATA11_1D)

   300  delete (DATA12_1D)

   301  delete (DATA22_1D)

   302 ;delete (range)

   303 ;delete (xvalues)

   304  delete (yvalues)

   305  delete (mn_yvalues)

   306  delete (mx_yvalues)

   307  delete (good)

   308  delete (max_bar)

   309  delete (min_bar)

   310  delete (max_cap)

   311  delete (min_cap)

   312 ;-----------------------------------------------------------------

   313 ;global res

   315   resg                     = True             ; Use plot options

   316   resg@cnFillOn            = True             ; Turn on color fill

   317   resg@gsnSpreadColors     = True             ; use full colormap

   318 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

   319 ; resg@lbLabelAutoStride   = True

   320   resg@cnLinesOn           = False            ; Turn off contourn lines

   321   resg@mpFillOn            = False            ; Turn off map fill

   323   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   324   resg@cnMinLevelValF       = 0.              ; Min level

   325   resg@cnMaxLevelValF       = 10.             ; Max level

   326   resg@cnLevelSpacingF      = 1.              ; interval

   328 ;global contour ob

   330   delta = 0.00001

   331   laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max)

   333   plot_name = "global_max_ob"

   334   title     = ob_name

   335   resg@tiMainString  = title

   337   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   338   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   340   plot = gsn_csm_contour_map_ce(wks,laiob_max,resg)

   341   frame(wks)

   343   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   344 ; system("rm "+plot_name+"."+plot_type)

   345 ; system("rm "+plot_name+"-1."+plot_type_new)

   346 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   348   clear (wks)

   349 ;------------------------------------------------------------------------

   350 ;global contour model

   352   plot_name = "global_max_model"

   353   title     = "Model " + model_name

   354   resg@tiMainString  = title

   356   wks = gsn_open_wks (plot_type,plot_name)

   357   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   359   delete (plot)

   360   plot = gsn_csm_contour_map_ce(wks,laimod_max,resg)

   361   frame(wks)

   363   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   364 ; system("rm "+plot_name+"."+plot_type)

   365 ; system("rm "+plot_name+"-1."+plot_type_new)

   366 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   368   clear (wks)

   369 ;------------------------------------------------------------------------

   370 ;global contour model vs ob

   372   plot_name = "global_max_model_vs_ob"

   374   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   375   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   377   delete (plot)

   378   plot=new(3,graphic)                        ; create graphic array

   380   resg@gsnFrame             = False          ; Do not draw plot

   381   resg@gsnDraw              = False          ; Do not advance frame

   383 ; plot correlation coef

   385   gRes               = True

   386   gRes@txFontHeightF = 0.02

   387   gRes@txAngleF      = 90

   389   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"

   391   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)

   392 ;--------------------------------------------------------------------

   394 ;(a) ob

   396   title     = ob_name

   397   resg@tiMainString  = title

   399   plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg)

   401 ;(b) model

   403   title     = "Model "+ model_name

   404   resg@tiMainString  = title

   406   plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg)

   408 ;(c) model-ob

   410   zz = laimod_max

   411   zz = laimod_max - laiob_max

   412   title = "Model_"+model_name+" - Observed"

   413   resg@tiMainString    = title

   415   resg@cnMinLevelValF  = -6.           ; Min level

   416   resg@cnMaxLevelValF  =  6.           ; Max level

   417   resg@cnLevelSpacingF =  1.           ; interval

   420   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)

   422   pres                            = True        ; panel plot mods desired

   423   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

   424                                                 ; indiv. plots in panel

   425   pres@gsnMaximize                = True        ; fill the page

   427   gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot

   429   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   430 ; system("rm "+plot_name+"."+plot_type)

   431 ; system("rm "+plot_name+"-1."+plot_type_new)

   432 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   434   frame (wks)

   435   clear (wks)

   437   delete (plot)

   438 ;-----------------------------------------------------------------

   439 ;(C) phase

   440 ;--------------------------------------------------------------------

   441 ; get data

   443 ; observed

   444   laiob_phase = laiob(0,:,:)

   445   s           = laiob(:,0,0)

   446   laiob_phase@long_name = "Leaf Area Index Max Month"

   448   dsizes_z = dimsizes(laiob)

   449   nlat     = dsizes_z(1)

   450   nlon     = dsizes_z(2)

   452   do j = 0,nlat-1

   453   do i = 0,nlon-1

   454      s = laiob(:,j,i)

   455      laiob_phase(j,i) = maxind(s) + 1

   456   end do

   457   end do

   459 ; print (min(laiob_phase)+"/"+max(laiob_phase))

   460   delete (s)

   461   delete (dsizes_z)

   462 ;-------------------------

   463 ; model

   464   laimod_phase = laimod(0,:,:)

   465   s            = laimod(:,0,0)

   466   laimod_phase@long_name = "Leaf Area Index Max Month"

   468   dsizes_z = dimsizes(laimod)

   469   nlat     = dsizes_z(1)

   470   nlon     = dsizes_z(2)

   472   do j = 0,nlat-1

   473   do i = 0,nlon-1

   474      s = laimod(:,j,i)

   475      laimod_phase(j,i) = maxind(s) + 1

   476   end do

   477   end do

   479 ; print (min(laimod_phase)+"/"+max(laimod_phase))

   480   delete (s)

   481   delete (dsizes_z)

   482 ;------------------------

   483   DATA11_1D = ndtooned(classob)

   484   DATA12_1D = ndtooned(laiob_phase)

   485   DATA22_1D = ndtooned(laimod_phase)

   487   yvalues      = new((/2,nx/),float)

   488   mn_yvalues   = new((/2,nx/),float)

   489   mx_yvalues   = new((/2,nx/),float)

   491   do nd=0,1

   493 ; See if we are doing model or observational data.

   495     if(nd.eq.0) then

   496       data_ob  = DATA11_1D

   497       data_mod = DATA12_1D

   498     else

   499       data_ob  = DATA11_1D

   500       data_mod = DATA22_1D

   501     end if

   503 ; Loop through each range and check for values.

   505     do i=0,nr-2

   506       if (i.ne.(nr-2)) then

   507 ;        print("")

   508 ;        print("In range ["+range(i)+","+range(i+1)+")")

   509         idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))

   510       else

   511 ;        print("")

   512 ;        print("In range ["+range(i)+",)")

   513         idx = ind(range(i).le.data_ob)

   514       end if

   516 ; Calculate average, and get min and max.

   518       if(.not.any(ismissing(idx))) then

   519         yvalues(nd,i)    = avg(data_mod(idx))

   520         mn_yvalues(nd,i) = min(data_mod(idx))

   521         mx_yvalues(nd,i) = max(data_mod(idx))

   522         count = dimsizes(idx)

   523       else

   524         count            = 0

   525         yvalues(nd,i)    = yvalues@_FillValue

   526         mn_yvalues(nd,i) = yvalues@_FillValue

   527         mx_yvalues(nd,i) = yvalues@_FillValue

   528       end if

   530 ;     print(nd + ": " + count + " points, avg = " + yvalues(nd,i))

   531 ;     print("Min/Max:  " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))

   533 ; Clean up for next time in loop.

   535       delete(idx)

   536     end do

   537     delete(data_ob)

   538     delete(data_mod)

   539   end do

   540 ;-----------------------------------------------------------------

   541 ; compute correlation coef and M score

   543   u = yvalues(0,:)

   544   v = yvalues(1,:)

   546   good = ind(.not.ismissing(u) .and. .not.ismissing(v))

   547   uu = u(good)

   548   vv = v(good)

   550   ccrPhase = esccr(uu,vv,0)

   551 ; print (ccrPhase)

   553 ; old eq

   554 ; bias   = abs(avg(vv)-avg(uu))

   555 ; new eq

   556   bias   = avg(abs(vv-uu))

   558   bias   = where((bias.gt. 6.),12.-bias,bias)

   559   Mphase = ((6. - bias)/6.)*5.

   561   print (Mphase)

   563  delete (u)

   564  delete (v)

   565  delete (uu)

   566  delete (vv)

   567 ;--------------------------------------------------------------------

   568 ; histogram res

   570   resm                = True

   571   resm@gsnMaximize    = True

   572   resm@gsnDraw        = False

   573   resm@gsnFrame       = False

   574   resm@xyMarkLineMode = "Markers"

   575   resm@xyMarkerSizeF  = 0.014

   576   resm@xyMarker       = 16

   577   resm@xyMarkerColors = (/"Brown","Blue"/)

   578 ; resm@trYMinF        = min(mn_yvalues) - 10.

   579 ; resm@trYMaxF        = max(mx_yvalues) + 10.

   580   resm@trYMinF        = min(mn_yvalues) - 2

   581   resm@trYMaxF        = max(mx_yvalues) + 4

   583   resm@tiYAxisString  = "Max LAI (Leaf Area Index) Month"

   584   resm@tiXAxisString  = "Land Cover Type"

   586   max_bar = new((/2,nx/),graphic)

   587   min_bar = new((/2,nx/),graphic)

   588   max_cap = new((/2,nx/),graphic)

   589   min_cap = new((/2,nx/),graphic)

   591   lnres = True

   592   line_colors = (/"brown","blue"/)

   593 ;------------------------------------------------------------------

   594 ; Start the graphics.

   596   plot_name = "histogram_phase"

   597   title     = model_name + " vs Observed"

   598   resm@tiMainString  = title

   600   wks   = gsn_open_wks (plot_type,plot_name)

   601 ;-----------------------------

   602 ; Add a boxed legend using the more simple method

   604   resm@pmLegendDisplayMode    = "Always"

   605 ; resm@pmLegendWidthF         = 0.1

   606   resm@pmLegendWidthF         = 0.08

   607   resm@pmLegendHeightF        = 0.05

   608   resm@pmLegendOrthogonalPosF = -1.17

   609 ; resm@pmLegendOrthogonalPosF = -1.00  ;(downward)

   610 ; resm@pmLegendParallelPosF   =  0.18

   611   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

   613 ; resm@lgPerimOn              = False

   614   resm@lgLabelFontHeightF     = 0.015

   615   resm@xyExplicitLegendLabels = (/"observed",model_name/)

   616 ;-----------------------------

   617   tRes  = True

   618   tRes@txFontHeightF = 0.025

   620   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"

   622   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

   624   xy = gsn_csm_xy(wks,xvalues,yvalues,resm)

   625 ;-------------------------------

   626 ;Attach the vertical bar and the horizontal cap line

   628   do nd=0,1

   629     lnres@gsLineColor = line_colors(nd)

   630     do i=0,nx-1

   632       if(.not.ismissing(mn_yvalues(nd,i)).and. \

   633          .not.ismissing(mx_yvalues(nd,i))) then

   635 ; Attach the vertical bar, both above and below the marker.

   637         x1 = xvalues(nd,i)

   638         y1 = yvalues(nd,i)

   639         y2 = mn_yvalues(nd,i)

   640         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

   642         y2 = mx_yvalues(nd,i)

   643         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

   645 ; Attach the horizontal cap line, both above and below the marker.

   647         x1 = xvalues(nd,i) - dx4

   648         x2 = xvalues(nd,i) + dx4

   649         y1 = mn_yvalues(nd,i)

   650         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

   652         y1 = mx_yvalues(nd,i)

   653         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

   654       end if

   655     end do

   656   end do

   658   draw(xy)

   659   frame(wks)

   661   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   662 ; system("rm "+plot_name+"."+plot_type)

   663 ; system("rm "+plot_name+"-1."+plot_type_new)

   664 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   666   clear (wks)

   668  delete (DATA11_1D)

   669  delete (DATA12_1D)

   670  delete (DATA22_1D)

   671 ;delete (range)

   672 ;delete (xvalues)

   673  delete (yvalues)

   674  delete (mn_yvalues)

   675  delete (mx_yvalues)

   676  delete (good)

   677  delete (max_bar)

   678  delete (min_bar)

   679  delete (max_cap)

   680  delete (min_cap)

   681 ;-----------------------------------------------------------------

   682 ;global res

   684   resg                     = True             ; Use plot options

   685   resg@cnFillOn            = True             ; Turn on color fill

   686   resg@gsnSpreadColors     = True             ; use full colormap

   687 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

   688 ; resg@lbLabelAutoStride   = True

   689   resg@cnLinesOn           = False            ; Turn off contourn lines

   690   resg@mpFillOn            = False            ; Turn off map fill

   692   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   693   resg@cnMinLevelValF       = 1.              ; Min level

   694   resg@cnMaxLevelValF       = 12.             ; Max level

   695   resg@cnLevelSpacingF      = 1.              ; interval

   697 ;global contour ob

   699   delta = 0.00001

   700   laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase)

   702   plot_name = "global_phase_ob"

   703   title     = ob_name

   704   resg@tiMainString  = title

   706   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   707   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   709   plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg)

   710   frame(wks)

   712   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   713 ; system("rm "+plot_name+"."+plot_type)

   714 ; system("rm "+plot_name+"-1."+plot_type_new)

   715 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   717   clear (wks)

   718 ;------------------------------------------------------------------------

   719 ;global contour model

   721   plot_name = "global_phase_model"

   722   title     = "Model " + model_name

   723   resg@tiMainString  = title

   725   wks = gsn_open_wks (plot_type,plot_name)

   726   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   728   delete (plot)

   729   plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg)

   730   frame(wks)

   732   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   733 ; system("rm "+plot_name+"."+plot_type)

   734 ; system("rm "+plot_name+"-1."+plot_type_new)

   735 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   737   clear (wks)

   738 ;------------------------------------------------------------------------

   739 ;global contour model vs ob

   741   plot_name = "global_phase_model_vs_ob"

   743   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   744   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   746   delete (plot)

   747   plot=new(3,graphic)                        ; create graphic array

   749   resg@gsnFrame             = False          ; Do not draw plot

   750   resg@gsnDraw              = False          ; Do not advance frame

   752 ; plot correlation coef

   754   gRes               = True

   755   gRes@txFontHeightF = 0.02

   756   gRes@txAngleF      = 90

   758   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"

   760   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)

   761 ;--------------------------------------------------------------------

   763 ;(a) ob

   765   title     = ob_name

   766   resg@tiMainString  = title

   768   plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg)

   770 ;(b) model

   772   title     = "Model "+ model_name

   773   resg@tiMainString  = title

   775   plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg)

   777 ;(c) model-ob

   779   delete (zz)

   780   zz = laimod_phase

   781   zz = laimod_phase - laiob_phase

   782   title = "Model_"+model_name+" - Observed"

   783   resg@tiMainString    = title

   785   resg@cnMinLevelValF  = -6.           ; Min level

   786   resg@cnMaxLevelValF  =  6.           ; Max level

   787   resg@cnLevelSpacingF =  1.           ; interval

   790   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)

   792 ; pres                            = True        ; panel plot mods desired

   793 ; pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

   794                                                 ; indiv. plots in panel

   795 ; pres@gsnMaximize                = True        ; fill the page

   797   gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot

   799   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   800 ; system("rm "+plot_name+"."+plot_type)

   801 ; system("rm "+plot_name+"-1."+plot_type_new)

   802 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   804   frame (wks)

   805   clear (wks)

   807   delete (plot)

   808 ;-----------------------------------------------------------------

   809 ;(D) grow day

   810 ;--------------------------------------------------------------------

   811 ; get data

   813   day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)

   815 ; observed

   816   laiob_grow = laiob(0,:,:)

   817   laiob_grow@long_name = "Days of Growing Season"

   819   dsizes_z = dimsizes(laiob)

   820   ntime    = dsizes_z(0)

   821   nlat     = dsizes_z(1)

   822   nlon     = dsizes_z(2)

   824   do j = 0,nlat-1

   825   do i = 0,nlon-1

   826      nday = 0.

   827      do k = 0,ntime-1

   828         if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then

   829            nday = nday + day_of_data(k)

   830         end if

   831      end do

   833      laiob_grow(j,i) = nday

   834   end do

   835   end do

   837 ; print (min(laiob_grow)+"/"+max(laiob_grow))

   838 ;-------------------------

   839 ; model

   840   laimod_grow = laimod(0,:,:)

   841   laimod_grow@long_name = "Days of Growing Season"

   843   dsizes_z = dimsizes(laimod)

   844   ntime    = dsizes_z(0)

   845   nlat     = dsizes_z(1)

   846   nlon     = dsizes_z(2)

   848   do j = 0,nlat-1

   849   do i = 0,nlon-1

   850      nday = 0.

   851      do k = 0,ntime-1

   852         if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then

   853            nday = nday + day_of_data(k)

   854         end if

   855      end do

   857      laimod_grow(j,i) = nday

   858   end do

   859   end do

   861 ; print (min(laimod_grow)+"/"+max(laimod_grow))

   862 ;------------------------

   863   DATA11_1D = ndtooned(classob)

   864   DATA12_1D = ndtooned(laiob_grow)

   865   DATA22_1D = ndtooned(laimod_grow)

   867   yvalues      = new((/2,nx/),float)

   868   mn_yvalues   = new((/2,nx/),float)

   869   mx_yvalues   = new((/2,nx/),float)

   871   do nd=0,1

   873 ; See if we are doing model or observational data.

   875     if(nd.eq.0) then

   876       data_ob  = DATA11_1D

   877       data_mod = DATA12_1D

   878     else

   879       data_ob  = DATA11_1D

   880       data_mod = DATA22_1D

   881     end if

   883 ; Loop through each range and check for values.

   885     do i=0,nr-2

   886       if (i.ne.(nr-2)) then

   887 ;        print("")

   888 ;        print("In range ["+range(i)+","+range(i+1)+")")

   889         idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))

   890       else

   891 ;        print("")

   892 ;        print("In range ["+range(i)+",)")

   893         idx = ind(range(i).le.data_ob)

   894       end if

   896 ; Calculate average, and get min and max.

   898       if(.not.any(ismissing(idx))) then

   899         yvalues(nd,i)    = avg(data_mod(idx))

   900         mn_yvalues(nd,i) = min(data_mod(idx))

   901         mx_yvalues(nd,i) = max(data_mod(idx))

   902         count = dimsizes(idx)

   903       else

   904         count            = 0

   905         yvalues(nd,i)    = yvalues@_FillValue

   906         mn_yvalues(nd,i) = yvalues@_FillValue

   907         mx_yvalues(nd,i) = yvalues@_FillValue

   908       end if

   910 ;     print(nd + ": " + count + " points, avg = " + yvalues(nd,i))

   911 ;     print("Min/Max:  " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))

   913 ; Clean up for next time in loop.

   915       delete(idx)

   916     end do

   917     delete(data_ob)

   918     delete(data_mod)

   919   end do

   920 ;-----------------------------------------------------------------

   921 ; compute correlation coef and M score

   923   u = yvalues(0,:)

   924   v = yvalues(1,:)

   926   good = ind(.not.ismissing(u) .and. .not.ismissing(v))

   927   uu = u(good)

   928   vv = v(good)

   930   ccrGrow = esccr(uu,vv,0)

   931 ; print (ccrGrow)

   933 ; new eq

   934   bias  = sum(abs(vv-uu)/(vv+uu))

   935   Mgrow = (1.- (bias/dimsizes(uu)))*5.

   937   print (Mgrow)

   939  delete (u)

   940  delete (v)

   941  delete (uu)

   942  delete (vv)

   943 ;--------------------------------------------------------------------

   944 ; histogram res

   946   resm                = True

   947   resm@gsnMaximize    = True

   948   resm@gsnDraw        = False

   949   resm@gsnFrame       = False

   950   resm@xyMarkLineMode = "Markers"

   951   resm@xyMarkerSizeF  = 0.014

   952   resm@xyMarker       = 16

   953   resm@xyMarkerColors = (/"Brown","Blue"/)

   954 ; resm@trYMinF        = min(mn_yvalues) - 10.

   955 ; resm@trYMaxF        = max(mx_yvalues) + 10.

   956   resm@trYMinF        = min(mn_yvalues) - 2

   957   resm@trYMaxF        = max(mx_yvalues) + 4

   959   resm@tiYAxisString = "Days of Growing season"

   960   resm@tiXAxisString = "Land Cover Type"

   962   max_bar = new((/2,nx/),graphic)

   963   min_bar = new((/2,nx/),graphic)

   964   max_cap = new((/2,nx/),graphic)

   965   min_cap = new((/2,nx/),graphic)

   967   lnres = True

   968   line_colors = (/"brown","blue"/)

   969 ;------------------------------------------------------------------

   970 ; Start the graphics.

   972   plot_name = "histogram_grow"

   973   title     = model_name + " vs Observed"

   974   resm@tiMainString  = title

   976   wks   = gsn_open_wks (plot_type,plot_name)

   977 ;-----------------------------

   978 ; Add a boxed legend using the more simple method

   980   resm@pmLegendDisplayMode    = "Always"

   981 ; resm@pmLegendWidthF         = 0.1

   982   resm@pmLegendWidthF         = 0.08

   983   resm@pmLegendHeightF        = 0.05

   984   resm@pmLegendOrthogonalPosF = -1.17

   985 ; resm@pmLegendOrthogonalPosF = -1.00  ;(downward)

   986 ; resm@pmLegendParallelPosF   =  0.18

   987   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

   989 ; resm@lgPerimOn              = False

   990   resm@lgLabelFontHeightF     = 0.015

   991   resm@xyExplicitLegendLabels = (/"observed",model_name/)

   992 ;-----------------------------

   993   tRes  = True

   994   tRes@txFontHeightF = 0.025

   996   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"

   998   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

  1000   xy = gsn_csm_xy(wks,xvalues,yvalues,resm)

  1001 ;-------------------------------

  1002 ;Attach the vertical bar and the horizontal cap line

  1004   do nd=0,1

  1005     lnres@gsLineColor = line_colors(nd)

  1006     do i=0,nx-1

  1008       if(.not.ismissing(mn_yvalues(nd,i)).and. \

  1009          .not.ismissing(mx_yvalues(nd,i))) then

  1011 ; Attach the vertical bar, both above and below the marker.

  1013         x1 = xvalues(nd,i)

  1014         y1 = yvalues(nd,i)

  1015         y2 = mn_yvalues(nd,i)

  1016         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1018         y2 = mx_yvalues(nd,i)

  1019         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1021 ; Attach the horizontal cap line, both above and below the marker.

  1023         x1 = xvalues(nd,i) - dx4

  1024         x2 = xvalues(nd,i) + dx4

  1025         y1 = mn_yvalues(nd,i)

  1026         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1028         y1 = mx_yvalues(nd,i)

  1029         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1030       end if

  1031     end do

  1032   end do

  1034   draw(xy)

  1035   frame(wks)

  1037   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1038 ; system("rm "+plot_name+"."+plot_type)

  1039 ; system("rm "+plot_name+"-1."+plot_type_new)

  1040 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1042   clear (wks)

  1044  delete (DATA11_1D)

  1045  delete (DATA12_1D)

  1046  delete (DATA22_1D)

  1047 ;delete (range)

  1048 ;delete (xvalues)

  1049  delete (yvalues)

  1050  delete (mn_yvalues)

  1051  delete (mx_yvalues)

  1052  delete (good)

  1053  delete (max_bar)

  1054  delete (min_bar)

  1055  delete (max_cap)

  1056  delete (min_cap)

  1057 ;-----------------------------------------------------------------

  1058 ;global res

  1060   resg                     = True             ; Use plot options

  1061   resg@cnFillOn            = True             ; Turn on color fill

  1062   resg@gsnSpreadColors     = True             ; use full colormap

  1063 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

  1064 ; resg@lbLabelAutoStride   = True

  1065   resg@cnLinesOn           = False            ; Turn off contourn lines

  1066   resg@mpFillOn            = False            ; Turn off map fill

  1068   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

  1069   resg@cnMinLevelValF       = 0.              ; Min level

  1070   resg@cnMaxLevelValF       = 360.            ; Max level

  1071   resg@cnLevelSpacingF      = 30.              ; interval

  1073 ;global contour ob

  1075   delta = 0.00001

  1076   laiob_grow = where(ismissing(laiob_grow).and.(ismissing(laimod_grow).or.(laimod_grow.lt.delta)),0.,laiob_grow)

  1078   plot_name = "global_grow_ob"

  1079   title     = ob_name

  1080   resg@tiMainString  = title

  1082   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1083   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1085   plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg)

  1086   frame(wks)

  1088   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1089 ; system("rm "+plot_name+"."+plot_type)

  1090 ; system("rm "+plot_name+"-1."+plot_type_new)

  1091 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1093   clear (wks)

  1094 ;------------------------------------------------------------------------

  1095 ;global contour model

  1097   plot_name = "global_grow_model"

  1098   title     = "Model " + model_name

  1099   resg@tiMainString  = title

  1101   wks = gsn_open_wks (plot_type,plot_name)

  1102   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1104   delete (plot)

  1105   plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg)

  1106   frame(wks)

  1108   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1109 ; system("rm "+plot_name+"."+plot_type)

  1110 ; system("rm "+plot_name+"-1."+plot_type_new)

  1111 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1113   clear (wks)

  1114 ;------------------------------------------------------------------------

  1115 ;global contour model vs ob

  1117   plot_name = "global_grow_model_vs_ob"

  1119   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1120   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1122   delete (plot)

  1123   plot=new(3,graphic)                        ; create graphic array

  1125   resg@gsnFrame             = False          ; Do not draw plot

  1126   resg@gsnDraw              = False          ; Do not advance frame

  1128 ; plot correlation coef

  1130   gRes               = True

  1131   gRes@txFontHeightF = 0.02

  1132   gRes@txAngleF      = 90

  1134   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"

  1136   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)

  1137 ;--------------------------------------------------------------------

  1139 ;(a) ob

  1141   title     = ob_name

  1142   resg@tiMainString  = title

  1144   plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg)

  1146 ;(b) model

  1148   title     = "Model "+ model_name

  1149   resg@tiMainString  = title

  1151   plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg)

  1153 ;(c) model-ob

  1155   delete (zz)

  1156   zz = laimod_grow

  1157   zz = laimod_grow - laiob_grow

  1158   title = "Model_"+model_name+" - Observed"

  1159   resg@tiMainString    = title

  1161   resg@cnMinLevelValF  = -120.           ; Min level

  1162   resg@cnMaxLevelValF  =  120.           ; Max level

  1163   resg@cnLevelSpacingF =  20.            ; interval

  1166   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)

  1168   pres                            = True        ; panel plot mods desired

  1169   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

  1170                                                 ; indiv. plots in panel

  1171   pres@gsnMaximize                = True        ; fill the page

  1173   gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot

  1175   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1176 ; system("rm "+plot_name+"."+plot_type)

  1177 ; system("rm "+plot_name+"-1."+plot_type_new)

  1178 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1180   frame (wks)

  1181   clear (wks)

  1183   delete (plot)

  1184 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 26 Jan 2009 22:08:20 -0500
changeset 0	0c6405ab2ff4
permissions	-rw-r--r--