clamp: npp/18.scatter_bias.ncl@0c6405ab2ff4

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

     2 ; combine  scatter, histogram, global and zonal plots

     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"

     7 ;************************************************

     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  plot_type     = "ps"

    20  plot_type_new = "png"

    22 ;************************************************

    23 ; read in ob data

    24 ;************************************************

    25 ;(1)

    26  dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

    27  fil81 = "data.81.nc"

    28  f81   = addfile (dir81+fil81,"r")

    30  id81   = f81->SITE_ID

    31  npp81  = f81->TNPP_C

    32  rain81 = tofloat(f81->PREC_ANN)

    33  x81    = f81->LONG_DD

    34  y81    = f81->LAT_DD

    36  id81@long_name  = "SITE_ID"

    38 ;change longitude from (-180,180) to (0,360)

    39 ;for model data interpolation

    41  do i= 0,dimsizes(x81)-1

    42     if (x81(i) .lt. 0.) then

    43         x81(i) = x81(i)+ 360.

    44     end if

    45  end do

    46 ;print (x81)

    47 ;-------------------------------------

    48 ;(2)

    49  dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

    50  fil933 = "data.933.nc"

    51  f933   = addfile (dir933+fil933,"r")

    53  id933   = f933->SITE_ID

    54  npp933  = f933->TNPP_C

    55  rain933 = f933->PREC

    56  x933    = f933->LONG_DD

    57  y933    = f933->LAT_DD

    59  id933@long_name  = "SITE_ID"

    61 ;change longitude from (-180,180) to (0,360)

    62 ;for model data interpolation

    64  do i= 0,dimsizes(x933)-1

    65     if (x933(i) .lt. 0.) then

    66         x933(i) = x933(i)+ 360.

    67     end if

    68  end do

    69 ;print (x933)

    70 ;----------------------------------------

    71 ;(3)

    72  dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

    73  filglobe = "Npp_T31_mean.nc"

    74  fglobe   = addfile (dirglobe+filglobe,"r")

    76  nppglobe0 = fglobe->NPP

    77  nppglobe  = nppglobe0

    78 ;************************************************

    79 ; read in model data

    80 ;************************************************

    81  model_name = "b30.061n"

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

    84  film = "b30.061n_1995-2004_ANN_climo_lnd.nc"

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

    87  nppmod0  = fm->NPP

    88  rainmod0 = fm->RAIN

    89  xm       = fm->lon

    90  ym       = fm->lat

    92  nppmod   = nppmod0(0,:,:)

    93  rainmod  = rainmod0(0,:,:)

    94  delete (nppmod0)

    95  delete (rainmod0)

    97  nppmod81  =linint2_points(xm,ym,nppmod,True,x81,y81,0)

    99  nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)

   101  rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)

   103  rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)

   105 ;************************************************

   106 ; Units for these variables are:

   107 ;

   108 ; rain81  : mm/year

   109 ; rainmod : mm/s

   110 ; npp81   : g C/m^2/year

   111 ; nppmod81: g C/m^2/s

   112 ; nppglobe: g C/m^2/year

   113 ;

   114 ; We want to convert these to "m/year" and "g C/m^2/year".

   115 ;

   116   nsec_per_year = 60*60*24*365

   118   rain81    = rain81 / 1000.

   119   rainmod81 = (rainmod81/ 1000.) * nsec_per_year

   120   nppmod81  = nppmod81 * nsec_per_year

   122   rain933    = rain933 / 1000.

   123   rainmod933 = (rainmod933/ 1000.) * nsec_per_year

   124   nppmod933  = nppmod933 * nsec_per_year

   126   nppmod  = nppmod * nsec_per_year

   128   npp81@units      = "gC/m^2/yr"

   129   nppmod81@units   = "gC/m^2/yr"

   130   npp933@units     = "gC/m^2/yr"

   131   nppmod933@units  = "gC/m^2/yr"

   132   nppmod@units     = "gC/m^2/yr"

   133   nppglobe@units   = "gC/m^2/yr"

   134   rain81@units     = "m/yr"

   135   rainmod81@units  = "m/yr"

   136   rain933@units    = "m/yr"

   137   rainmod933@units = "m/yr"

   139   npp81@long_name      = "NPP (gC/m2/year)"

   140   npp933@long_name     = "NPP (gC/m2/year)"

   141   nppmod81@long_name   = "NPP (gC/m2/year)"

   142   nppmod933@long_name  = "NPP (gC/m2/year)"

   143   nppmod@long_name     = "NPP (gC/m2/year)"

   144   nppglobe@long_name   = "NPP (gC/m2/year)"

   145   rain81@long_name     = "PREC (m/year)"

   146   rain933@long_name    = "PREC (m/year)"

   147   rainmod81@long_name  = "PREC (m/year)"

   148   rainmod933@long_name = "PREC (m/year)"

   150 ;(A) plot scatter ob 81

   152  plot_name = "scatter_ob_81"

   153  title     = "Observed 81 sites"

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

   156  res                   = True                  ; plot mods desired

   157  res@tiMainString      = title                 ; add title

   158  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   159  res@xyMarkers         =  16                   ; choose type of marker

   160  res@xyMarkerColor     = "red"                 ; Marker color

   161  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   162  res@tmLabelAutoStride = True                  ; nice tick mark labels

   164  plot  = gsn_csm_xy (wks,id81,npp81,res)       ; create plot

   165  frame(wks)

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

   168  system("rm "+plot_name+"."+plot_type)

   169  system("rm "+plot_name+"-1."+plot_type_new)

   170  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   172  clear (wks)

   174 ;(B) plot scatter ob 933

   176  plot_name = "scatter_ob_933"

   177  title     = "Observed 933 sites"

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

   180  res                   = True                  ; plot mods desired

   181  res@tiMainString      = title                 ; add title

   182  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   183  res@xyMarkers         =  16                   ; choose type of marker

   184  res@xyMarkerColor     = "red"                 ; Marker color

   185  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   186  res@tmLabelAutoStride = True                  ; nice tick mark labels

   188  plot  = gsn_csm_xy (wks,id933,npp933,res)     ; create plot

   189  frame(wks)

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

   192  system("rm "+plot_name+"."+plot_type)

   193  system("rm "+plot_name+"-1."+plot_type_new)

   194  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   196  clear (wks)

   198 ;(C) plot scatter model 81

   200  plot_name = "scatter_model_81"

   201  title     = "Model "+ model_name +" 81 sites"

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

   204  res                   = True                  ; plot mods desired

   205  res@tiMainString      = title                 ; add title

   206  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   207  res@xyMarkers         =  16                   ; choose type of marker

   208  res@xyMarkerColor     = "red"                 ; Marker color

   209  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   210  res@tmLabelAutoStride = True                  ; nice tick mark labels

   212  plot  = gsn_csm_xy (wks,id81,nppmod81,res)    ; create plot

   213  frame(wks)

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

   216  system("rm "+plot_name+"."+plot_type)

   217  system("rm "+plot_name+"-1."+plot_type_new)

   218  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   220  clear (wks)

   222 ;(D) plot scatter model 933

   224  plot_name = "scatter_model_933"

   225  title     = "Model "+ model_name +" 933 sites"

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

   228  res                   = True                  ; plot mods desired

   229  res@tiMainString      = title                 ; add title

   230  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   231  res@xyMarkers         =  16                   ; choose type of marker

   232  res@xyMarkerColor     = "red"                 ; Marker color

   233  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   234  res@tmLabelAutoStride = True                  ; nice tick mark labels

   236  plot  = gsn_csm_xy (wks,id933,nppmod933,res)    ; create plot

   237  frame(wks)

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

   240  system("rm "+plot_name+"."+plot_type)

   241  system("rm "+plot_name+"-1."+plot_type_new)

   242  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   244  clear (wks)

   246 ;(E) scatter model vs ob 81

   248  plot_name = "scatter_model_vs_ob_81"

   249  title     = model_name +" vs ob 81 sites"

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

   252  res                   = True                  ; plot mods desired

   253  res@tiMainString      = title                 ; add title

   254  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   255  res@xyMarkers         =  16                   ; choose type of marker

   256  res@xyMarkerColor     = "red"                 ; Marker color

   257  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   258  res@tmLabelAutoStride = True                  ; nice tick mark labels

   259 ;-------------------------------

   260 ;compute correlation coef. and M

   261  ccr81 = esccr(nppmod81,npp81,0)

   262  print (ccr81)

   264 ;new eq

   265  bias = sum(abs(nppmod81-npp81)/(nppmod81+npp81))

   266  M81  = (1. - (bias/dimsizes(y81)))*5.

   267  print (M81)

   269  tRes  = True

   270  tRes@txFontHeightF = 0.025

   272  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"

   274  gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)

   275 ;--------------------------------

   276  plot  = gsn_csm_xy (wks,npp81,nppmod81,res)       ; create plot

   277  frame(wks)

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

   280  system("rm "+plot_name+"."+plot_type)

   281  system("rm "+plot_name+"-1."+plot_type_new)

   282  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   284  clear (wks)

   286 ;(F) scatter model vs ob 933

   288  plot_name = "scatter_model_vs_ob_933"

   289  title     = model_name +" vs ob 933 sites"

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

   292  res                   = True                  ; plot mods desired

   293  res@tiMainString      = title                 ; add title

   294  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   295  res@xyMarkers         =  16                   ; choose type of marker

   296  res@xyMarkerColor     = "red"                 ; Marker color

   297  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   298  res@tmLabelAutoStride = True                  ; nice tick mark labels

   299 ;-------------------------------

   300 ;compute correlation coef. and M

   301  ccr933 = esccr(nppmod933,npp933,0)

   302  print (ccr933)

   304 ;new eq

   305  bias = sum(abs(nppmod933-npp933)/(nppmod933+npp933))

   306  M933 = (1. - (bias/dimsizes(y933)))*5.

   307  print (M933)

   309  tRes  = True

   310  tRes@txFontHeightF = 0.025

   312  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"

   314  gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)

   315 ;--------------------------------

   316  plot  = gsn_csm_xy (wks,npp933,nppmod933,res)    ; create plot

   317  frame(wks)

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

   320  system("rm "+plot_name+"."+plot_type)

   321  system("rm "+plot_name+"-1."+plot_type_new)

   322  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   324  clear (wks)

   325 ;**************************************************************************

   326 ;(G) histogram 81

   328 ;--------------------------------------------------------------------

   329 ;get data

   331   RAIN1_1D = ndtooned(rain81)

   332   RAIN2_1D = ndtooned(rainmod81)

   333   NPP1_1D  = ndtooned(npp81)

   334   NPP2_1D  = ndtooned(nppmod81)

   336 ; Calculaee "nice" bins for binning the data in equally spaced ranges.

   338   nbins       = 15     ; Number of bins to use.

   340   nicevals    = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)

   341   nvals       = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)

   342   range       = fspan(nicevals(0),nicevals(1),nvals)

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

   345 ; particular range, and then take an average.

   347   nr           = dimsizes(range)

   348   nx           = nr-1

   349   xvalues      = new((/2,nx/),typeof(RAIN1_1D))

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

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

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

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

   354   yvalues      = new((/2,nx/),typeof(RAIN1_1D))

   355   mn_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

   356   mx_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

   358   do nd=0,1

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

   362     if(nd.eq.0) then

   363       data     = RAIN1_1D

   364       npp_data = NPP1_1D

   365     else

   366       data     = RAIN2_1D

   367       npp_data = NPP2_1D

   368     end if

   370 ; Loop through each range and check for values.

   372     do i=0,nr-2

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

   374 ;        print("")

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

   376         idx = ind((range(i).le.data).and.(data.lt.range(i+1)))

   377       else

   378 ;        print("")

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

   380         idx = ind(range(i).le.data)

   381       end if

   383 ; Calculate average, and get min and max.

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

   386         yvalues(nd,i)    = avg(npp_data(idx))

   387         mn_yvalues(nd,i) = min(npp_data(idx))

   388         mx_yvalues(nd,i) = max(npp_data(idx))

   389         count = dimsizes(idx)

   390       else

   391         count            = 0

   392         yvalues(nd,i)    = yvalues@_FillValue

   393         mn_yvalues(nd,i) = yvalues@_FillValue

   394         mx_yvalues(nd,i) = yvalues@_FillValue

   395       end if

   397 ; Print out information.

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

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

   403 ; Clean up for next time in loop.

   405       delete(idx)

   406     end do

   407     delete(data)

   408     delete(npp_data)

   409   end do

   410 ;----------------------------------------

   411 ;compute correlation coeff and M score

   413  u = yvalues(0,:)

   414  v = yvalues(1,:)

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

   417  uu = u(good)

   418  vv = v(good)

   420  ccr81h = esccr(uu,vv,0)

   421  print (ccr81h)

   423 ;new eq

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

   425  M81h = (1.- (bias/dimsizes(uu)))*5.

   426  print (M81h)

   428  delete (u)

   429  delete (v)

   430  delete (uu)

   431  delete (vv)

   432 ;----------------------------------------------------------------------

   433 ; histogram res

   435   resh                = True

   436   resh@gsnMaximize    = True

   437   resh@gsnDraw        = False

   438   resh@gsnFrame       = False

   439   resh@xyMarkLineMode = "Markers"

   440   resh@xyMarkerSizeF  = 0.014

   441   resh@xyMarker       = 16

   442   resh@xyMarkerColors = (/"Brown","Blue"/)

   443   resh@trYMinF        = min(mn_yvalues) - 10.

   444   resh@trYMaxF        = max(mx_yvalues) + 10.

   446   resh@tiYAxisString  = "NPP (g C/m2/year)"

   447   resh@tiXAxisString  = "Precipitation (m/year)"

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

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

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

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

   454   lnres = True

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

   456 ;=================================================================

   457 ; histogram ob 81 site only

   458 ;

   459   plot_name = "histogram_ob_81"

   460   title     = "Observed 81 site"

   461   resh@tiMainString  = title

   463   wks   = gsn_open_wks (plot_type,plot_name)

   465   xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)

   467 ;-------------------------------

   468 ;Attach the vertical bar and the horizontal cap line

   470   do nd=0,0

   471     lnres@gsLineColor = line_colors(nd)

   472     do i=0,nx-1

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

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

   476 ;

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

   478 ;

   479         x1 = xvalues(nd,i)

   480         y1 = yvalues(nd,i)

   481         y2 = mn_yvalues(nd,i)

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

   484         y2 = mx_yvalues(nd,i)

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

   486 ;

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

   488 ;

   489         x1 = xvalues(nd,i) - dx4

   490         x2 = xvalues(nd,i) + dx4

   491         y1 = mn_yvalues(nd,i)

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

   494         y1 = mx_yvalues(nd,i)

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

   496       end if

   497     end do

   498   end do

   500   draw(xy)

   501   frame(wks)

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

   504   system("rm "+plot_name+"."+plot_type)

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

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

   508   clear (wks)

   509 ;===========================================================================

   510 ; histogram model vs ob 81 site

   512   plot_name = "histogram_mod-ob_81"

   513   title     = model_name+ " vs Observed 81 site"

   514   resh@tiMainString  = title

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

   518 ;-----------------------------

   519 ; Add a boxed legend using the more simple method, which won't have

   520 ; vertical lines going through the markers.

   522   resh@pmLegendDisplayMode    = "Always"

   523 ; resh@pmLegendWidthF         = 0.1

   524   resh@pmLegendWidthF         = 0.08

   525   resh@pmLegendHeightF        = 0.05

   526   resh@pmLegendOrthogonalPosF = -1.17

   527 ; resh@pmLegendOrthogonalPosF = -1.00  ;(downward)

   528 ; resh@pmLegendParallelPosF   =  0.18

   529   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

   531 ; resh@lgPerimOn              = False

   532   resh@lgLabelFontHeightF     = 0.015

   533   resh@xyExplicitLegendLabels = (/"observed",model_name/)

   534 ;-----------------------------

   535   tRes  = True

   536   tRes@txFontHeightF = 0.025

   538   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"

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

   542   xy = gsn_csm_xy(wks,xvalues,yvalues,resh)

   543 ;-------------------------------

   544 ;Attach the vertical bar and the horizontal cap line

   546   do nd=0,1

   547     lnres@gsLineColor = line_colors(nd)

   548     do i=0,nx-1

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

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

   552 ;

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

   554 ;

   555         x1 = xvalues(nd,i)

   556         y1 = yvalues(nd,i)

   557         y2 = mn_yvalues(nd,i)

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

   560         y2 = mx_yvalues(nd,i)

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

   562 ;

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

   564 ;

   565         x1 = xvalues(nd,i) - dx4

   566         x2 = xvalues(nd,i) + dx4

   567         y1 = mn_yvalues(nd,i)

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

   570         y1 = mx_yvalues(nd,i)

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

   572       end if

   573     end do

   574   end do

   575   draw(xy)

   576   frame(wks)

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

   579   system("rm "+plot_name+"."+plot_type)

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

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

   583  clear (wks)

   585  delete (RAIN1_1D)

   586  delete (RAIN2_1D)

   587  delete (NPP1_1D)

   588  delete (NPP2_1D)

   589  delete (range)

   590  delete (xvalues)

   591  delete (yvalues)

   592  delete (mn_yvalues)

   593  delete (mx_yvalues)

   594  delete (good)

   595  delete (max_bar)

   596  delete (min_bar)

   597  delete (max_cap)

   598  delete (min_cap)

   599 ;**************************************************************************

   600 ;(H) histogram 933

   602 ;--------------------------------------------------------------------

   603 ;get data

   605   RAIN1_1D = ndtooned(rain933)

   606   RAIN2_1D = ndtooned(rainmod933)

   607   NPP1_1D  = ndtooned(npp933)

   608   NPP2_1D  = ndtooned(nppmod933)

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

   612   nbins       = 15     ; Number of bins to use.

   614   nicevals    = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)

   615   nvals       = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)

   616   range       = fspan(nicevals(0),nicevals(1),nvals)

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

   619 ; particular range, and then take an average.

   621   nr           = dimsizes(range)

   622   nx           = nr-1

   623   xvalues      = new((/2,nx/),typeof(RAIN1_1D))

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

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

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

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

   628   yvalues      = new((/2,nx/),typeof(RAIN1_1D))

   629   mn_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

   630   mx_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

   632   do nd=0,1

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

   636     if(nd.eq.0) then

   637       data     = RAIN1_1D

   638       npp_data = NPP1_1D

   639     else

   640       data     = RAIN2_1D

   641       npp_data = NPP2_1D

   642     end if

   644 ; Loop through each range and check for values.

   646     do i=0,nr-2

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

   648 ;        print("")

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

   650         idx = ind((range(i).le.data).and.(data.lt.range(i+1)))

   651       else

   652 ;        print("")

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

   654         idx = ind(range(i).le.data)

   655       end if

   657 ; Calculate average, and get min and max.

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

   660         yvalues(nd,i)    = avg(npp_data(idx))

   661         mn_yvalues(nd,i) = min(npp_data(idx))

   662         mx_yvalues(nd,i) = max(npp_data(idx))

   663         count = dimsizes(idx)

   664       else

   665         count            = 0

   666         yvalues(nd,i)    = yvalues@_FillValue

   667         mn_yvalues(nd,i) = yvalues@_FillValue

   668         mx_yvalues(nd,i) = yvalues@_FillValue

   669       end if

   671 ; Print out information.

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

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

   676 ; Clean up for next time in loop.

   678       delete(idx)

   679     end do

   680     delete(data)

   681     delete(npp_data)

   682   end do

   683 ;----------------------------------------

   684 ;compute correlation coeff and M score

   686  u = yvalues(0,:)

   687  v = yvalues(1,:)

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

   690  uu = u(good)

   691  vv = v(good)

   693  ccr933h = esccr(uu,vv,0)

   694  print (ccr933h)

   696 ;new eq

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

   698  M933h = (1.- (bias/dimsizes(uu)))*5.

   699  print (M933h)

   701  delete (u)

   702  delete (v)

   703  delete (uu)

   704  delete (vv)

   705 ;----------------------------------------------------------------------

   706 ; histogram res

   708   delete (resh)

   709   resh                = True

   710   resh@gsnMaximize    = True

   711   resh@gsnDraw        = False

   712   resh@gsnFrame       = False

   713   resh@xyMarkLineMode = "Markers"

   714   resh@xyMarkerSizeF  = 0.014

   715   resh@xyMarker       = 16

   716   resh@xyMarkerColors = (/"Brown","Blue"/)

   717   resh@trYMinF        = min(mn_yvalues) - 10.

   718   resh@trYMaxF        = max(mx_yvalues) + 10.

   720   resh@tiYAxisString  = "NPP (g C/m2/year)"

   721   resh@tiXAxisString  = "Precipitation (m/year)"

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

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

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

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

   728   lnres = True

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

   730 ;=================================================================

   731 ; histogram ob 933 site only

   733   plot_name = "histogram_ob_933"

   734   title     = "Observed 933 site"

   735   resh@tiMainString  = title

   737   wks   = gsn_open_wks (plot_type,plot_name)

   739   xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)

   741 ;-------------------------------

   742 ;Attach the vertical bar and the horizontal cap line

   744   do nd=0,0

   745     lnres@gsLineColor = line_colors(nd)

   746     do i=0,nx-1

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

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

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

   753         x1 = xvalues(nd,i)

   754         y1 = yvalues(nd,i)

   755         y2 = mn_yvalues(nd,i)

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

   758         y2 = mx_yvalues(nd,i)

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

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

   763         x1 = xvalues(nd,i) - dx4

   764         x2 = xvalues(nd,i) + dx4

   765         y1 = mn_yvalues(nd,i)

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

   768         y1 = mx_yvalues(nd,i)

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

   770       end if

   771     end do

   772   end do

   774   draw(xy)

   775   frame(wks)

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

   778   system("rm "+plot_name+"."+plot_type)

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

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

   782   delete (xy)

   783   clear (wks)

   785 ;===========================================================================

   786 ; histogram model vs ob 933 site

   788   plot_name = "histogram_mod-ob_933"

   789   title     = model_name+ " vs Observed 933 site"

   790   resh@tiMainString  = title

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

   794 ;-----------------------------

   795 ; Add a boxed legend using the more simple method, which won't have

   796 ; vertical lines going through the markers.

   798   resh@pmLegendDisplayMode    = "Always"

   799 ; resh@pmLegendWidthF         = 0.1

   800   resh@pmLegendWidthF         = 0.08

   801   resh@pmLegendHeightF        = 0.05

   802   resh@pmLegendOrthogonalPosF = -1.17

   803 ; resh@pmLegendOrthogonalPosF = -1.00  ;(downward)

   804 ; resh@pmLegendParallelPosF   =  0.18

   805   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

   807 ; resh@lgPerimOn              = False

   808   resh@lgLabelFontHeightF     = 0.015

   809   resh@xyExplicitLegendLabels = (/"observed",model_name/)

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

   811   tRes  = True

   812   tRes@txFontHeightF = 0.025

   814   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"

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

   818   xy = gsn_csm_xy(wks,xvalues,yvalues,resh)

   819 ;-------------------------------

   820 ;Attach the vertical bar and the horizontal cap line

   822   do nd=0,1

   823     lnres@gsLineColor = line_colors(nd)

   824     do i=0,nx-1

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

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

   828 ;

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

   830 ;

   831         x1 = xvalues(nd,i)

   832         y1 = yvalues(nd,i)

   833         y2 = mn_yvalues(nd,i)

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

   836         y2 = mx_yvalues(nd,i)

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

   838 ;

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

   840 ;

   841         x1 = xvalues(nd,i) - dx4

   842         x2 = xvalues(nd,i) + dx4

   843         y1 = mn_yvalues(nd,i)

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

   846         y1 = mx_yvalues(nd,i)

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

   848       end if

   849     end do

   850   end do

   851   draw(xy)

   852   frame(wks)

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

   855   system("rm "+plot_name+"."+plot_type)

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

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

   859   delete (xy)

   860   clear (wks)

   861 ;------------------------------------------------------------------------

   862 ; global contour

   864 ;res

   865   resg                     = True             ; Use plot options

   866   resg@cnFillOn            = True             ; Turn on color fill

   867   resg@gsnSpreadColors     = True             ; use full colormap

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

   869 ; resg@lbLabelAutoStride   = True

   870   resg@cnLinesOn           = False            ; Turn off contourn lines

   871   resg@mpFillOn            = False            ; Turn off map fill

   873   resg@gsnSpreadColors      = True            ; use full colormap

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

   875   resg@cnMinLevelValF       = 0.              ; Min level

   876   resg@cnMaxLevelValF       = 2200.           ; Max level

   877   resg@cnLevelSpacingF      = 200.            ; interval

   878 ;------------------------------------------------------------------------

   879 ;global contour ob

   881   delta = 0.00000000001

   882   nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)

   884   plot_name = "global_ob"

   885   title     = "Observed MODIS MOD 17"

   886   resg@tiMainString  = title

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

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

   891   plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)

   892   frame(wks)

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

   895   system("rm "+plot_name+"."+plot_type)

   896   system("rm "+plot_name+"-1."+plot_type_new)

   897   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   899   clear (wks)

   900 ;------------------------------------------------------------------------

   901 ;global contour model

   903   plot_name = "global_model"

   904   title     = "Model "+ model_name

   905   resg@tiMainString  = title

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

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

   910   plot = gsn_csm_contour_map_ce(wks,nppmod,resg)

   911   frame(wks)

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

   914   system("rm "+plot_name+"."+plot_type)

   915   system("rm "+plot_name+"-1."+plot_type_new)

   916   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   918   clear (wks)

   919 ;------------------------------------------------------------------------

   920 ;global contour model vs ob

   922   plot_name = "global_model_vs_ob"

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

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

   927   delete (plot)

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

   930   resg@gsnFrame             = False          ; Do not draw plot

   931   resg@gsnDraw              = False          ; Do not advance frame

   933 ;(d) compute correlation coef and M score

   935   uu1 = ndtooned(nppmod)

   936   vv1 = ndtooned(nppglobe)

   938   delete (good)

   939   good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))

   941   ug = uu1(good)

   942   vg = vv1(good)

   944   ccrG = esccr(ug,vg,0)

   945   print (ccrG)

   947   MG   = (ccrG*ccrG)* 5.0

   948   print (MG)

   950 ; plot correlation coef

   952   gRes  = True

   953   gRes@txFontHeightF = 0.02

   954   gRes@txAngleF = 90

   956   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"

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

   959 ;--------------------------------------------------------------------

   961 ;(a) ob

   963   title     = "Observed MODIS MOD 17"

   964   resg@tiMainString  = title

   966   plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)

   968 ;(b) model

   970   title     = "Model "+ model_name

   971   resg@tiMainString  = title

   973   plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)

   975 ;(c) model-ob

   977   zz = nppmod

   978   zz = nppmod - nppglobe

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

   981   resg@cnMinLevelValF  = -500           ; Min level

   982   resg@cnMaxLevelValF  =  500.          ; Max level

   983   resg@cnLevelSpacingF =  50.           ; interval

   984   resg@tiMainString    = title

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

   988   pres                            = True        ; panel plot mods desired

   989   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

   990                                                 ; indiv. plots in panel

   991   pres@gsnMaximize                = True        ; fill the page

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

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

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

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

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

  1000   frame (wks)

  1001   clear (wks)

  1003   delete (plot)

  1004 ;---------------------------------------------------------------------

  1005 ; zonal line plot ob

  1007   resz = True

  1009   vv     = zonalAve(nppglobe)

  1010   vv@long_name = nppglobe@long_name

  1012   plot_name = "zonal_ob"

  1013   title     = "Observed MODIS MOD 17"

  1014   resz@tiMainString  = title

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

  1018   plot  = gsn_csm_xy (wks,ym,vv,resz)

  1019   frame(wks)

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

  1022   system("rm "+plot_name+"."+plot_type)

  1023   system("rm "+plot_name+"-1."+plot_type_new)

  1024   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1026   clear (wks)

  1027 ;---------------------------------------------------------------------

  1028 ; zonal line plot model vs ob

  1030   s  = new ((/2,dimsizes(ym)/), float)

  1032   s(0,:) = zonalAve(nppglobe)

  1033   s(1,:) = zonalAve(nppmod)

  1035   s@long_name = nppglobe@long_name

  1036 ;-------------------------------------------

  1037 ;(d) compute correlation coef and M score

  1039   ccrZ = esccr(s(1,:), s(0,:),0)

  1040   print (ccrZ)

  1042   MZ   = (ccrZ*ccrZ)* 5.0

  1043   print (MZ)

  1044 ;-------------------------------------------

  1045   plot_name = "zonal_model_vs_ob"

  1046   title     = "Zonal Average"

  1047   resz@tiMainString  = title

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

  1051 ; resz@vpHeightF          = 0.4               ; change aspect ratio of plot

  1052 ; resz@vpWidthF           = 0.7

  1054   resz@xyMonoLineColor    = "False"           ; want colored lines

  1055   resz@xyLineColors       = (/"black","red"/) ; colors chosen

  1056 ; resz@xyLineThicknesses  = (/3.,3./)      ; line thicknesses

  1057   resz@xyLineThicknesses  = (/2.,2./)      ; line thicknesses

  1058   resz@xyDashPatterns     = (/0.,0./)      ; make all lines solid

  1060   resz@tiYAxisString      = s@long_name      ; add a axis title

  1061   resz@txFontHeightF      = 0.0195            ; change title font heights

  1063 ; Legent

  1064   resz@pmLegendDisplayMode    = "Always"            ; turn on legend

  1065   resz@pmLegendSide           = "Top"               ; Change location of

  1066 ; resz@pmLegendParallelPosF   = .45                 ; move units right

  1067   resz@pmLegendParallelPosF   = .82                 ; move units right

  1068   resz@pmLegendOrthogonalPosF = -0.4                ; move units down

  1070   resz@pmLegendWidthF         = 0.10                ; Change width and

  1071   resz@pmLegendHeightF        = 0.10                ; height of legend.

  1072   resz@lgLabelFontHeightF     = .02                 ; change font height

  1073 ; resz@lgTitleOn              = True                ; turn on legend title

  1074 ; resz@lgTitleString          = "Example"           ; create legend title

  1075 ; resz@lgTitleFontHeightF     = .025                ; font of legend title

  1076   resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels

  1077 ;--------------------------------------------------------------------

  1078   zRes  = True

  1079   zRes@txFontHeightF = 0.025

  1081   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"

  1083   gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)

  1084 ;--------------------------------------------------------------------

  1086   plot  = gsn_csm_xy (wks,ym,s,resz)       ; create plot

  1088   frame(wks)                                            ; advance frame

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

  1091   system("rm "+plot_name+"."+plot_type)

  1092   system("rm "+plot_name+"-1."+plot_type_new)

  1093   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1095   clear (wks)

  1096 end

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