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

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

     2 ; combine all scatter

     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  dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

    26  fil81 = "data.81.nc"

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

    29  id81   = f81->SITE_ID

    30  npp81  = f81->TNPP_C

    31  rain81 = tofloat(f81->PREC_ANN)

    32  x81    = f81->LONG_DD

    33  y81    = f81->LAT_DD

    35  id81@long_name  = "SITE_ID"

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

    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  dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

    49  fil933 = "data.933.nc"

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

    52  id933   = f933->SITE_ID

    53  npp933  = f933->TNPP_C

    54  rain933 = f933->PREC

    55  x933    = f933->LONG_DD

    56  y933    = f933->LAT_DD

    58  id933@long_name  = "SITE_ID"

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

    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 ; read in model data

    72 ;************************************************

    73  model_name = "b30.061n"

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

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

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

    79  nppmod  = fm->NPP

    80  rainmod = fm->RAIN

    81  xm      = fm->lon

    82  ym      = fm->lat

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

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

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

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

    92 ;************************************************

    93 ; Units for these four variables are:

    94 ;

    95 ; rain81  : mm/year

    96 ; rainmod : mm/s

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

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

    99 ;

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

   101 ;

   102   nsec_per_year = 60*60*24*365

   104   rain81    = rain81 / 1000.

   105   rainmod81 = (rainmod81/ 1000.) * nsec_per_year

   106   nppmod81  = nppmod81 * nsec_per_year

   108   rain933    = rain933 / 1000.

   109   rainmod933 = (rainmod933/ 1000.) * nsec_per_year

   110   nppmod933  = nppmod933 * nsec_per_year

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

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

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

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

   116   rain81@units     = "m/yr"

   117   rainmod81@units  = "m/yr"

   118   rain933@units    = "m/yr"

   119   rainmod933@units = "m/yr"

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

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

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

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

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

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

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

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

   130 ;(A) plot scatter ob 81

   132  plot_name = "scatter_ob_81"

   133  title     = "Observed 81 sites"

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

   136  res                   = True                  ; plot mods desired

   137  res@tiMainString      = title                 ; add title

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

   139  res@xyMarkers         =  16                   ; choose type of marker

   140  res@xyMarkerColor     = "red"                 ; Marker color

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

   142  res@tmLabelAutoStride = True                  ; nice tick mark labels

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

   145  frame(wks)

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

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

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

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

   152  clear (wks)

   154 ;(B) plot scatter ob 933

   156  plot_name = "scatter_ob_933"

   157  title     = "Observed 933 sites"

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

   160  res                   = True                  ; plot mods desired

   161  res@tiMainString      = title                 ; add title

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

   163  res@xyMarkers         =  16                   ; choose type of marker

   164  res@xyMarkerColor     = "red"                 ; Marker color

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

   166  res@tmLabelAutoStride = True                  ; nice tick mark labels

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

   169  frame(wks)

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

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

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

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

   176  clear (wks)

   178 ;(C) plot scatter model 81

   180  plot_name = "scatter_model_81"

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

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

   184  res                   = True                  ; plot mods desired

   185  res@tiMainString      = title                 ; add title

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

   187  res@xyMarkers         =  16                   ; choose type of marker

   188  res@xyMarkerColor     = "red"                 ; Marker color

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

   190  res@tmLabelAutoStride = True                  ; nice tick mark labels

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

   193  frame(wks)

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

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

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

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

   200  clear (wks)

   202 ;(D) plot scatter model 933

   204  plot_name = "scatter_model_933"

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

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

   208  res                   = True                  ; plot mods desired

   209  res@tiMainString      = title                 ; add title

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

   211  res@xyMarkers         =  16                   ; choose type of marker

   212  res@xyMarkerColor     = "red"                 ; Marker color

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

   214  res@tmLabelAutoStride = True                  ; nice tick mark labels

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

   217  frame(wks)

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

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

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

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

   224  clear (wks)

   226 ;(E) scatter model vs ob 81

   228  plot_name = "scatter_model_vs_ob_81"

   229  title     = model_name +" vs ob 81 sites"

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

   232  res                   = True                  ; plot mods desired

   233  res@tiMainString      = title                 ; add title

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

   235  res@xyMarkers         =  16                   ; choose type of marker

   236  res@xyMarkerColor     = "red"                 ; Marker color

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

   238  res@tmLabelAutoStride = True                  ; nice tick mark labels

   239 ;-------------------------------

   240 ;compute correlation coef. and M

   241  ccr81 = esccr(nppmod81,npp81,0)

   242  print (ccr81)

   244 ;new eq

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

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

   247  print (M81)

   249  tRes  = True

   250  tRes@txFontHeightF = 0.025

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

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

   255 ;--------------------------------

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

   257  frame(wks)

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

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

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

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

   264  clear (wks)

   266 ;(F) scatter model vs ob 933

   268  plot_name = "scatter_model_vs_ob_933"

   269  title     = model_name +" vs ob 933 sites"

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

   272  res                   = True                  ; plot mods desired

   273  res@tiMainString      = title                 ; add title

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

   275  res@xyMarkers         =  16                   ; choose type of marker

   276  res@xyMarkerColor     = "red"                 ; Marker color

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

   278  res@tmLabelAutoStride = True                  ; nice tick mark labels

   279 ;-------------------------------

   280 ;compute correlation coef. and M

   281  ccr933 = esccr(nppmod933,npp933,0)

   282  print (ccr933)

   284 ;new eq

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

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

   287  print (M933)

   289  tRes  = True

   290  tRes@txFontHeightF = 0.025

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

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

   295 ;--------------------------------

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

   297  frame(wks)

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

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

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

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

   304  clear (wks)

   305 ;**************************************************************************

   306 ;(G) histogram 81

   308 ;--------------------------------------------------------------------

   309 ;get data

   311   RAIN1_1D = ndtooned(rain81)

   312   RAIN2_1D = ndtooned(rainmod81)

   313   NPP1_1D  = ndtooned(npp81)

   314   NPP2_1D  = ndtooned(nppmod81)

   315 ;

   316 ; Calculate some "nice" bins for binning the data in equally spaced

   317 ; ranges.

   318 ;

   319   nbins       = 15     ; Number of bins to use.

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

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

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

   324 ;

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

   326 ; particular range, and then take an average.

   327 ;

   328   nr           = dimsizes(range)

   329   nx           = nr-1

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

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

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

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

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

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

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

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

   339   do nd=0,1

   340 ;

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

   342 ;

   343     if(nd.eq.0) then

   344       data     = RAIN1_1D

   345       npp_data = NPP1_1D

   346     else

   347       data     = RAIN2_1D

   348       npp_data = NPP2_1D

   349     end if

   350 ;

   351 ; Loop through each range and check for values.

   352 ;

   353     do i=0,nr-2

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

   355 ;        print("")

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

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

   358       else

   359 ;        print("")

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

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

   362       end if

   363 ;

   364 ; Calculate average, and get min and max.

   365 ;

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

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

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

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

   370         count = dimsizes(idx)

   371       else

   372         count            = 0

   373         yvalues(nd,i)    = yvalues@_FillValue

   374         mn_yvalues(nd,i) = yvalues@_FillValue

   375         mx_yvalues(nd,i) = yvalues@_FillValue

   376       end if

   377 ;

   378 ; Print out information.

   379 ;

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

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

   383 ;

   384 ; Clean up for next time in loop.

   385 ;

   386       delete(idx)

   387     end do

   388     delete(data)

   389     delete(npp_data)

   390   end do

   391 ;----------------------------------------

   392 ;compute correlation coeff and M score

   393  u = yvalues(0,:)

   394  v = yvalues(1,:)

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

   397  uu = u(good)

   398  vv = v(good)

   400  ccr81h = esccr(uu,vv,0)

   401  print (ccr81h)

   403 ;new eq

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

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

   406  print (M81h)

   408  delete (u)

   409  delete (v)

   410  delete (uu)

   411  delete (vv)

   412 ;----------------------------------------------------------------------

   413 ; histogram res

   415   res                = True

   416   res@gsnMaximize    = True

   417   res@gsnDraw        = False

   418   res@gsnFrame       = False

   419   res@xyMarkLineMode = "Markers"

   420   res@xyMarkerSizeF  = 0.014

   421   res@xyMarker       = 16

   422   res@xyMarkerColors = (/"Brown","Blue"/)

   423   res@trYMinF        = min(mn_yvalues) - 10.

   424   res@trYMaxF        = max(mx_yvalues) + 10.

   426   res@tiYAxisString  = "NPP (g C/m2/year)"

   427   res@tiXAxisString  = "Precipitation (m/year)"

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

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

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

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

   434   lnres = True

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

   436 ;=================================================================

   437 ; histogram ob 81 site only

   438 ;

   439   plot_name = "histogram_ob_81"

   440   title     = "Observed 81 site"

   441   res@tiMainString  = title

   443   wks   = gsn_open_wks (plot_type,plot_name)

   445   xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),res)

   447 ;-------------------------------

   448 ;Attach the vertical bar and the horizontal cap line

   450   do nd=0,0

   451     lnres@gsLineColor = line_colors(nd)

   452     do i=0,nx-1

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

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

   456 ;

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

   458 ;

   459         x1 = xvalues(nd,i)

   460         y1 = yvalues(nd,i)

   461         y2 = mn_yvalues(nd,i)

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

   464         y2 = mx_yvalues(nd,i)

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

   466 ;

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

   468 ;

   469         x1 = xvalues(nd,i) - dx4

   470         x2 = xvalues(nd,i) + dx4

   471         y1 = mn_yvalues(nd,i)

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

   474         y1 = mx_yvalues(nd,i)

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

   476       end if

   477     end do

   478   end do

   480   draw(xy)

   481   frame(wks)

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

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

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

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

   488   clear (wks)

   489 ;===========================================================================

   490 ; histogram model vs ob 81 site

   492   plot_name = "histogram_mod-ob_81"

   493   title     = model_name+ " vs Observed 81 site"

   494   res@tiMainString  = title

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

   498 ;-----------------------------

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

   500 ; vertical lines going through the markers.

   502   res@pmLegendDisplayMode    = "Always"

   503 ; res@pmLegendWidthF         = 0.1

   504   res@pmLegendWidthF         = 0.08

   505   res@pmLegendHeightF        = 0.05

   506   res@pmLegendOrthogonalPosF = -1.17

   507 ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

   508 ; res@pmLegendParallelPosF   =  0.18

   509   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   511 ; res@lgPerimOn              = False

   512   res@lgLabelFontHeightF     = 0.015

   513   res@xyExplicitLegendLabels = (/"observed",model_name/)

   514 ;-----------------------------

   515   tRes  = True

   516   tRes@txFontHeightF = 0.025

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

   520   gsn_text_ndc(wks,correlation_text,0.5,0.8,tRes)

   522   xy = gsn_csm_xy(wks,xvalues,yvalues,res)

   523 ;-------------------------------

   524 ;Attach the vertical bar and the horizontal cap line

   526   do nd=0,1

   527     lnres@gsLineColor = line_colors(nd)

   528     do i=0,nx-1

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

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

   532 ;

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

   534 ;

   535         x1 = xvalues(nd,i)

   536         y1 = yvalues(nd,i)

   537         y2 = mn_yvalues(nd,i)

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

   540         y2 = mx_yvalues(nd,i)

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

   542 ;

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

   544 ;

   545         x1 = xvalues(nd,i) - dx4

   546         x2 = xvalues(nd,i) + dx4

   547         y1 = mn_yvalues(nd,i)

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

   550         y1 = mx_yvalues(nd,i)

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

   552       end if

   553     end do

   554   end do

   555   draw(xy)

   556   frame(wks)

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

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

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

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

   563  clear (wks)

   565  delete (RAIN1_1D)

   566  delete (RAIN2_1D)

   567  delete (NPP1_1D)

   568  delete (NPP2_1D)

   569  delete (range)

   570  delete (xvalues)

   571  delete (yvalues)

   572  delete (mn_yvalues)

   573  delete (mx_yvalues)

   574  delete (good)

   575  delete (max_bar)

   576  delete (min_bar)

   577  delete (max_cap)

   578  delete (min_cap)

   579 ;**************************************************************************

   580 ;(H) histogram 933

   582 ;--------------------------------------------------------------------

   583 ;get data

   585   RAIN1_1D = ndtooned(rain933)

   586   RAIN2_1D = ndtooned(rainmod933)

   587   NPP1_1D  = ndtooned(npp933)

   588   NPP2_1D  = ndtooned(nppmod933)

   589 ;

   590 ; Calculate some "nice" bins for binning the data in equally spaced

   591 ; ranges.

   592 ;

   593   nbins       = 15     ; Number of bins to use.

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

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

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

   598 ;

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

   600 ; particular range, and then take an average.

   601 ;

   602   nr           = dimsizes(range)

   603   nx           = nr-1

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

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

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

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

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

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

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

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

   613   do nd=0,1

   614 ;

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

   616 ;

   617     if(nd.eq.0) then

   618       data     = RAIN1_1D

   619       npp_data = NPP1_1D

   620     else

   621       data     = RAIN2_1D

   622       npp_data = NPP2_1D

   623     end if

   624 ;

   625 ; Loop through each range and check for values.

   626 ;

   627     do i=0,nr-2

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

   629 ;        print("")

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

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

   632       else

   633 ;        print("")

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

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

   636       end if

   637 ;

   638 ; Calculate average, and get min and max.

   639 ;

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

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

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

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

   644         count = dimsizes(idx)

   645       else

   646         count            = 0

   647         yvalues(nd,i)    = yvalues@_FillValue

   648         mn_yvalues(nd,i) = yvalues@_FillValue

   649         mx_yvalues(nd,i) = yvalues@_FillValue

   650       end if

   651 ;

   652 ; Print out information.

   653 ;

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

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

   657 ;

   658 ; Clean up for next time in loop.

   659 ;

   660       delete(idx)

   661     end do

   662     delete(data)

   663     delete(npp_data)

   664   end do

   665 ;----------------------------------------

   666 ;compute correlation coeff and M score

   667  u = yvalues(0,:)

   668  v = yvalues(1,:)

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

   671  uu = u(good)

   672  vv = v(good)

   674  ccr933h = esccr(uu,vv,0)

   675  print (ccr933h)

   677 ;new eq

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

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

   680  print (M933h)

   682  delete (u)

   683  delete (v)

   684  delete (uu)

   685  delete (vv)

   686 ;----------------------------------------------------------------------

   687 ; histogram res

   689   res                = True

   690   res@gsnMaximize    = True

   691   res@gsnDraw        = False

   692   res@gsnFrame       = False

   693   res@xyMarkLineMode = "Markers"

   694   res@xyMarkerSizeF  = 0.014

   695   res@xyMarker       = 16

   696   res@xyMarkerColors = (/"Brown","Blue"/)

   697   res@trYMinF        = min(mn_yvalues) - 10.

   698   res@trYMaxF        = max(mx_yvalues) + 10.

   700   res@tiYAxisString  = "NPP (g C/m2/year)"

   701   res@tiXAxisString  = "Precipitation (m/year)"

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

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

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

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

   708   lnres = True

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

   710 ;=================================================================

   711 ; histogram ob 933 site only

   712 ;

   713   plot_name = "histogram_ob_933"

   714   title     = "Observed 933 site"

   715   res@tiMainString  = title

   717   wks   = gsn_open_wks (plot_type,plot_name)

   719   xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),res)

   721 ;-------------------------------

   722 ;Attach the vertical bar and the horizontal cap line

   724   do nd=0,0

   725     lnres@gsLineColor = line_colors(nd)

   726     do i=0,nx-1

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

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

   730 ;

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

   732 ;

   733         x1 = xvalues(nd,i)

   734         y1 = yvalues(nd,i)

   735         y2 = mn_yvalues(nd,i)

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

   738         y2 = mx_yvalues(nd,i)

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

   740 ;

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

   742 ;

   743         x1 = xvalues(nd,i) - dx4

   744         x2 = xvalues(nd,i) + dx4

   745         y1 = mn_yvalues(nd,i)

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

   748         y1 = mx_yvalues(nd,i)

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

   750       end if

   751     end do

   752   end do

   754   draw(xy)

   755   frame(wks)

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

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

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

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

   762   clear (wks)

   763 ;===========================================================================

   764 ; histogram model vs ob 933 site

   766   plot_name = "histogram_mod-ob_933"

   767   title     = model_name+ " vs Observed 933 site"

   768   res@tiMainString  = title

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

   772 ;-----------------------------

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

   774 ; vertical lines going through the markers.

   776   res@pmLegendDisplayMode    = "Always"

   777 ; res@pmLegendWidthF         = 0.1

   778   res@pmLegendWidthF         = 0.08

   779   res@pmLegendHeightF        = 0.05

   780   res@pmLegendOrthogonalPosF = -1.17

   781 ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

   782 ; res@pmLegendParallelPosF   =  0.18

   783   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   785 ; res@lgPerimOn              = False

   786   res@lgLabelFontHeightF     = 0.015

   787   res@xyExplicitLegendLabels = (/"observed",model_name/)

   788 ;-----------------------------

   789   tRes  = True

   790   tRes@txFontHeightF = 0.025

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

   794   gsn_text_ndc(wks,correlation_text,0.5,0.8,tRes)

   796   xy = gsn_csm_xy(wks,xvalues,yvalues,res)

   797 ;-------------------------------

   798 ;Attach the vertical bar and the horizontal cap line

   800   do nd=0,1

   801     lnres@gsLineColor = line_colors(nd)

   802     do i=0,nx-1

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

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

   806 ;

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

   808 ;

   809         x1 = xvalues(nd,i)

   810         y1 = yvalues(nd,i)

   811         y2 = mn_yvalues(nd,i)

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

   814         y2 = mx_yvalues(nd,i)

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

   816 ;

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

   818 ;

   819         x1 = xvalues(nd,i) - dx4

   820         x2 = xvalues(nd,i) + dx4

   821         y1 = mn_yvalues(nd,i)

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

   824         y1 = mx_yvalues(nd,i)

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

   826       end if

   827     end do

   828   end do

   829   draw(xy)

   830   frame(wks)

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

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

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

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

   837  clear (wks)

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

   840 end

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