clamp: npp/17.scatter_bias.ncl@4be95183fbcd

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

     2 ; combine all scatter and all histogram

     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  nppglobe  = fglobe->NPP

    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)

   335 ;

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

   337 ; ranges.

   338 ;

   339   nbins       = 15     ; Number of bins to use.

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

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

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

   344 ;

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

   346 ; particular range, and then take an average.

   347 ;

   348   nr           = dimsizes(range)

   349   nx           = nr-1

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

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

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

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

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

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

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

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

   359   do nd=0,1

   360 ;

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

   362 ;

   363     if(nd.eq.0) then

   364       data     = RAIN1_1D

   365       npp_data = NPP1_1D

   366     else

   367       data     = RAIN2_1D

   368       npp_data = NPP2_1D

   369     end if

   370 ;

   371 ; Loop through each range and check for values.

   372 ;

   373     do i=0,nr-2

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

   375 ;        print("")

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

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

   378       else

   379 ;        print("")

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

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

   382       end if

   383 ;

   384 ; Calculate average, and get min and max.

   385 ;

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

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

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

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

   390         count = dimsizes(idx)

   391       else

   392         count            = 0

   393         yvalues(nd,i)    = yvalues@_FillValue

   394         mn_yvalues(nd,i) = yvalues@_FillValue

   395         mx_yvalues(nd,i) = yvalues@_FillValue

   396       end if

   397 ;

   398 ; Print out information.

   399 ;

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

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

   403 ;

   404 ; Clean up for next time in loop.

   405 ;

   406       delete(idx)

   407     end do

   408     delete(data)

   409     delete(npp_data)

   410   end do

   411 ;----------------------------------------

   412 ;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   res                = True

   436   res@gsnMaximize    = True

   437   res@gsnDraw        = False

   438   res@gsnFrame       = False

   439   res@xyMarkLineMode = "Markers"

   440   res@xyMarkerSizeF  = 0.014

   441   res@xyMarker       = 16

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

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

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

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

   447   res@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   res@tiMainString  = title

   463   wks   = gsn_open_wks (plot_type,plot_name)

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

   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   res@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   res@pmLegendDisplayMode    = "Always"

   523 ; res@pmLegendWidthF         = 0.1

   524   res@pmLegendWidthF         = 0.08

   525   res@pmLegendHeightF        = 0.05

   526   res@pmLegendOrthogonalPosF = -1.17

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

   528 ; res@pmLegendParallelPosF   =  0.18

   529   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   531 ; res@lgPerimOn              = False

   532   res@lgLabelFontHeightF     = 0.015

   533   res@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.5,0.8,tRes)

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

   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)

   609 ;

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

   611 ; ranges.

   612 ;

   613   nbins       = 15     ; Number of bins to use.

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

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

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

   618 ;

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

   620 ; particular range, and then take an average.

   621 ;

   622   nr           = dimsizes(range)

   623   nx           = nr-1

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

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

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

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

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

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

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

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

   633   do nd=0,1

   634 ;

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

   636 ;

   637     if(nd.eq.0) then

   638       data     = RAIN1_1D

   639       npp_data = NPP1_1D

   640     else

   641       data     = RAIN2_1D

   642       npp_data = NPP2_1D

   643     end if

   644 ;

   645 ; Loop through each range and check for values.

   646 ;

   647     do i=0,nr-2

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

   649 ;        print("")

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

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

   652       else

   653 ;        print("")

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

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

   656       end if

   657 ;

   658 ; Calculate average, and get min and max.

   659 ;

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

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

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

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

   664         count = dimsizes(idx)

   665       else

   666         count            = 0

   667         yvalues(nd,i)    = yvalues@_FillValue

   668         mn_yvalues(nd,i) = yvalues@_FillValue

   669         mx_yvalues(nd,i) = yvalues@_FillValue

   670       end if

   671 ;

   672 ; Print out information.

   673 ;

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

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

   677 ;

   678 ; Clean up for next time in loop.

   679 ;

   680       delete(idx)

   681     end do

   682     delete(data)

   683     delete(npp_data)

   684   end do

   685 ;----------------------------------------

   686 ;compute correlation coeff and M score

   687  u = yvalues(0,:)

   688  v = yvalues(1,:)

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

   691  uu = u(good)

   692  vv = v(good)

   694  ccr933h = esccr(uu,vv,0)

   695  print (ccr933h)

   697 ;new eq

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

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

   700  print (M933h)

   702  delete (u)

   703  delete (v)

   704  delete (uu)

   705  delete (vv)

   706 ;----------------------------------------------------------------------

   707 ; histogram res

   709   res                = True

   710   res@gsnMaximize    = True

   711   res@gsnDraw        = False

   712   res@gsnFrame       = False

   713   res@xyMarkLineMode = "Markers"

   714   res@xyMarkerSizeF  = 0.014

   715   res@xyMarker       = 16

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

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

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

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

   721   res@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

   732 ;

   733   plot_name = "histogram_ob_933"

   734   title     = "Observed 933 site"

   735   res@tiMainString  = title

   737   wks   = gsn_open_wks (plot_type,plot_name)

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

   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

   750 ;

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

   752 ;

   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)

   760 ;

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

   762 ;

   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   clear (wks)

   783 ;===========================================================================

   784 ; histogram model vs ob 933 site

   786   plot_name = "histogram_mod-ob_933"

   787   title     = model_name+ " vs Observed 933 site"

   788   res@tiMainString  = title

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

   792 ;-----------------------------

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

   794 ; vertical lines going through the markers.

   796   res@pmLegendDisplayMode    = "Always"

   797 ; res@pmLegendWidthF         = 0.1

   798   res@pmLegendWidthF         = 0.08

   799   res@pmLegendHeightF        = 0.05

   800   res@pmLegendOrthogonalPosF = -1.17

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

   802 ; res@pmLegendParallelPosF   =  0.18

   803   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   805 ; res@lgPerimOn              = False

   806   res@lgLabelFontHeightF     = 0.015

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

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

   809   tRes  = True

   810   tRes@txFontHeightF = 0.025

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

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

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

   817 ;-------------------------------

   818 ;Attach the vertical bar and the horizontal cap line

   820   do nd=0,1

   821     lnres@gsLineColor = line_colors(nd)

   822     do i=0,nx-1

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

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

   826 ;

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

   828 ;

   829         x1 = xvalues(nd,i)

   830         y1 = yvalues(nd,i)

   831         y2 = mn_yvalues(nd,i)

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

   834         y2 = mx_yvalues(nd,i)

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

   836 ;

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

   838 ;

   839         x1 = xvalues(nd,i) - dx4

   840         x2 = xvalues(nd,i) + dx4

   841         y1 = mn_yvalues(nd,i)

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

   844         y1 = mx_yvalues(nd,i)

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

   846       end if

   847     end do

   848   end do

   849   draw(xy)

   850   frame(wks)

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

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

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

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

   857   clear (wks)

   858 ;------------------------------------------------------------------------

   859 ;global contour

   861 ;res

   862   resw                     = True             ; Use plot options

   863   resw@cnFillOn            = True             ; Turn on color fill

   864   resw@gsnSpreadColors     = True             ; use full colormap

   865 ; resw@cnFillMode          = "RasterFill"     ; Turn on raster color

   866 ; resw@lbLabelAutoStride   = True

   867   resw@cnLinesOn           = False            ; Turn off contourn lines

   868   resw@mpFillOn            = False            ; Turn off map fill

   870   resw@gsnSpreadColors      = True            ; use full colormap

   871   resw@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   872   resw@cnMinLevelValF       = 0.              ; Min level

   873   resw@cnMaxLevelValF       = 2200.           ; Max level

   874   resw@cnLevelSpacingF      = 200.            ; interval

   875 ;------------------------------------------------------------------------

   876 ;global contour ob

   878   delta = 0.00000000001

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

   881   plot_name = "global_ob"

   882   title     = "Observed MODIS MOD 17"

   883   resw@tiMainString  = title

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

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

   889   plot = gsn_csm_contour_map_ce(wks,nppglobe,resw)

   890   frame(wks)

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

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

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

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

   897   clear (wks)

   898 ;------------------------------------------------------------------------

   899 ;global contour model

   901   plot_name = "global_model"

   902   title     = "Model "+ model_name

   903   resw@tiMainString  = title

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

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

   908   plot = gsn_csm_contour_map_ce(wks,nppmod,resw)

   909   frame(wks)

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

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

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

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

   916   clear (wks)

   917 ;------------------------------------------------------------------------

   918 ;global contour model vs ob

   920   plot_name = "global_model_vs_ob"

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

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

   925   delete (plot)

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

   928   resw@gsnFrame             = False          ; Do not draw plot

   929   resw@gsnDraw              = False          ; Do not advance frame

   931 ;(a) ob

   933   title     = "Observed MODIS MOD 17"

   934   resw@tiMainString  = title

   936   plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resw)        ; for observed

   938 ;(b) model

   940   title     = "Model "+ model_name

   941   resw@tiMainString  = title

   943   plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resw) ; for model

   945 ;(c) model-ob

   947   zz = nppmod

   948   zz = nppmod - nppglobe

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

   951   resw@cnMinLevelValF  = -500           ; Min level

   952   resw@cnMaxLevelValF  =  500.          ; Max level

   953   resw@cnLevelSpacingF =  50.           ; interval

   954   resw@tiMainString    = title

   956   plot(2) = gsn_csm_contour_map_ce(wks,zz,resw)

   958   pres                            = True        ; panel plot mods desired

   959   pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

   960                                                 ; indiv. plots in panel

   961   pres@gsnMaximize                = True        ; fill the page

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

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

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

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

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

   970   frame (wks)

   971   clear (wks)

   972 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--