clamp: npp/99.all.ncl.1@0c6405ab2ff4

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

     2 ; combine scatter, histogram, global and zonal plots

     3 ; compute all correlation coef and M score

     4 ; add 1-to-1 line in scatter plots

     5 ;*****************************************************

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

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

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

     9 ;*****************************************************

    10 procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \

    11                  ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \

    12                  ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \

    13                  ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \

    14                  ,dx4[1]:numeric \

    15                   )

    16 begin

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

    18 ; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D

    19 ; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4

    21   nbins       = 15     ; Number of bins to use.

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

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

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

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

    28 ; particular range, and then take an average.

    30   nr           = dimsizes(range)

    31   nx           = nr-1

    32 ; print (nx)

    34 ; xvalues      = new((/2,nx/),typeof(RAIN1_1D))

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

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

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

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

    39 ; yvalues      = new((/2,nx/),typeof(RAIN1_1D))

    40 ; mn_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

    41 ; mx_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

    43   do nd=0,1

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

    47     if(nd.eq.0) then

    48       data     = RAIN1_1D

    49       npp_data = NPP1_1D

    50     else

    51       data     = RAIN2_1D

    52       npp_data = NPP2_1D

    53     end if

    55 ; Loop through each range and check for values.

    57     do i=0,nr-2

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

    59 ;        print("")

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

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

    62       else

    63 ;        print("")

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

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

    66       end if

    68 ; Calculate average, and get min and max.

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

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

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

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

    74         count = dimsizes(idx)

    75       else

    76         count            = 0

    77         yvalues(nd,i)    = yvalues@_FillValue

    78         mn_yvalues(nd,i) = yvalues@_FillValue

    79         mx_yvalues(nd,i) = yvalues@_FillValue

    80       end if

    82 ; Print out information.

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

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

    88 ; Clean up for next time in loop.

    90       delete(idx)

    91     end do

    92     delete(data)

    93     delete(npp_data)

    94   end do

    95 end

    96 ;****************************************************************************

    97 ; Main code.

    98 ;****************************************************************************

    99 begin

   101  plot_type     = "ps"

   102  plot_type_new = "png"

   104 ;************************************************

   105 ; read model data

   106 ;************************************************

   108 ;film = "i01.06cn_1798-2004_ANN_climo.nc"

   109 ;model_name = "06cn"

   110 ;model_grid = "T42"

   112 ;film = "i01.06casa_1798-2004_ANN_climo.nc"

   113 ;model_name = "06casa"

   114 ;model_grid = "T42"

   116  film = "i01.10casa_1948-2004_ANN_climo.nc"

   117  model_name = "10casa"

   118  model_grid = "T42"

   120 ;film = "i01.10cn_1948-2004_ANN_climo.nc"

   121 ;model_name = "10cn"

   122 ;model_grid = "T42"

   123 ;--------------------------------------------------

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

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

   127  nppmod0  = fm->NPP

   128  rainmod0 = fm->RAIN

   129  xm       = fm->lon

   130  ym       = fm->lat

   132 ;************************************************

   133 ; read ob data

   134 ;************************************************

   136 ;(1) data at 81 sites

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

   138  fil81 = "data.81.nc"

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

   141  id81   = f81->SITE_ID

   142  npp81  = f81->TNPP_C

   143  rain81 = tofloat(f81->PREC_ANN)

   144  x81    = f81->LONG_DD

   145  y81    = f81->LAT_DD

   147  id81@long_name  = "SITE_ID"

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

   150 ;for model data interpolation

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

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

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

   155     end if

   156  end do

   157 ;print (x81)

   158 ;-------------------------------------

   159 ;(2) data at 933 sites

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

   161  fil933 = "data.933.nc"

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

   164  id933   = f933->SITE_ID

   165  npp933  = f933->TNPP_C

   166  rain933 = f933->PREC

   167  x933    = f933->LONG_DD

   168  y933    = f933->LAT_DD

   170  id933@long_name  = "SITE_ID"

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

   173 ;for model data interpolation

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

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

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

   178     end if

   179  end do

   180 ;print (x933)

   181 ;----------------------------------------

   182 ;(3) global data, interpolated into model grid

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

   184  filglobe = "Npp_"+model_grid+"_mean.nc"

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

   187  nppglobe0 = fglobe->NPP

   188  nppglobe  = nppglobe0

   190 ;***********************************************************************

   191 ; interpolate model data into ob sites

   192 ;***********************************************************************

   194  nppmod   = nppmod0(0,:,:)

   195  rainmod  = rainmod0(0,:,:)

   196  delete (nppmod0)

   197  delete (rainmod0)

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

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

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

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

   207 ;**********************************************************

   208 ; unified units

   209 ;**********************************************************

   210 ; Units for these variables are:

   211 ;

   212 ; rain81  : mm/year

   213 ; rainmod : mm/s

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

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

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

   217 ;

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

   220   nsec_per_year = 60*60*24*365

   222   rain81    = rain81 / 1000.

   223   rainmod81 = (rainmod81/ 1000.) * nsec_per_year

   224   nppmod81  = nppmod81 * nsec_per_year

   226   rain933    = rain933 / 1000.

   227   rainmod933 = (rainmod933/ 1000.) * nsec_per_year

   228   nppmod933  = nppmod933 * nsec_per_year

   230   nppmod  = nppmod * nsec_per_year

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

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

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

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

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

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

   238   rain81@units     = "m/yr"

   239   rainmod81@units  = "m/yr"

   240   rain933@units    = "m/yr"

   241   rainmod933@units = "m/yr"

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

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

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

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

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

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

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

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

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

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

   254 ;*******************************************************************

   255 ;(A)-1 table of site81 (1) (the first 41 sites)

   256 ;*******************************************************************

   258   table_length = 0.95

   260 ; table header name

   261   table_header_name = "Site ID"

   263 ; column (not including header column)

   264   col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)

   265   ncol = dimsizes(col_header)

   267 ; row (not including header row)

   268   nrow       = 41

   269   row_header = new ((/nrow/),string )

   270   row_header(0:nrow-1) = id81(0:nrow-1)

   272 ; Table header

   273   ncr1  = (/1,1/)               ; 1 row, 1 column

   274   x1    = (/0.005,0.15/)        ; Start and end X

   275   y1    = (/0.900,0.95/)       ; Start and end Y

   276   text1 = table_header_name

   277   res1               = True

   278   res1@txFontHeightF = 0.015

   279   res1@gsFillColor   = "CornFlowerBlue"

   281 ; Column header (equally space in x2)

   282   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   283   x2    = (/x1(1),0.995/)       ; start from end of x1

   284   y2    = y1                    ; same as y1

   285   text2 = col_header

   286   res2               = True

   287   res2@txFontHeightF = 0.015

   288   res2@gsFillColor   = "Gray"

   290 ; Row header (equally space in y2)

   291   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   292   x3    = x1                         ; same as x1

   293   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   294   text3 = row_header

   295   res3               = True

   296   res3@txFontHeightF = 0.010

   297   res3@gsFillColor   = "Gray"

   299 ; Main table body

   300   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   301   x4    = x2                      ; Start and end x

   302   y4    = y3                      ; Start and end Y

   303   text4 = new((/nrow,ncol/),string)

   305   color_fill4           = new((/nrow,ncol/),string)

   306   color_fill4           = "white"

   307 ; color_fill4(:,ncol-1) = "grey"

   308 ; color_fill4(nrow-1,:) = "green"

   310   res4               = True       ; Set up resource list

   311 ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

   312   res4@txFontHeightF = 0.015

   313   res4@gsFillColor   = color_fill4

   315   delete (color_fill4)

   316 ;-------------------------------------------------------------------

   317 ; for table value

   319   do n = 0,nrow-1

   320      text4(n,0) = sprintf("%5.2f", y81(n))

   321      text4(n,1) = sprintf("%5.2f", x81(n))

   322      text4(n,2) = sprintf("%5.2f", npp81(n))

   323      text4(n,3) = sprintf("%5.2f", rain81(n))

   324   end do

   325 ;-------------------------------------------------------------------

   327   plot_name = "table_site81_ob1"

   329   wks = gsn_open_wks (plot_type,plot_name)

   330 ;------------------------------------------

   331 ; for table title

   333   gRes               = True

   334   gRes@txFontHeightF = 0.02

   335 ; gRes@txAngleF      = 90

   337   title_text = "Observation at 81 Sites (1)"

   339   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   340 ;------------------------------------------

   342   gsn_table(wks,ncr1,x1,y1,text1,res1)

   343   gsn_table(wks,ncr2,x2,y2,text2,res2)

   344   gsn_table(wks,ncr3,x3,y3,text3,res3)

   345   gsn_table(wks,ncr4,x4,y4,text4,res4)

   347   frame(wks)

   348   clear (wks)

   350   delete (row_header)

   351   delete (res3)

   352   delete (ncr3)

   353   delete (text3)

   354   delete (res4)

   355   delete (ncr4)

   356   delete (text4)

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

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

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

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

   363 ;*******************************************************************

   364 ;(A)-2 table of site81 (2) (the last 40 sites)

   365 ;*******************************************************************

   367   table_length = 0.95

   369 ; table header name

   370   table_header_name = "Site ID"

   372 ; column (not including header column)

   373   col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)

   374   ncol = dimsizes(col_header)

   376 ; row (not including header row)

   377   nrow       = 40

   378   row_header = new ((/nrow/),string )

   379   row_header(0:nrow-1) = id81(41:41+nrow-1)

   381 ; Table header

   382   ncr1  = (/1,1/)               ; 1 row, 1 column

   383   x1    = (/0.005,0.15/)        ; Start and end X

   384   y1    = (/0.900,0.95/)       ; Start and end Y

   385   text1 = table_header_name

   386   res1               = True

   387   res1@txFontHeightF = 0.015

   388   res1@gsFillColor   = "CornFlowerBlue"

   390 ; Column header (equally space in x2)

   391   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   392   x2    = (/x1(1),0.995/)       ; start from end of x1

   393   y2    = y1                    ; same as y1

   394   text2 = col_header

   395   res2               = True

   396   res2@txFontHeightF = 0.015

   397   res2@gsFillColor   = "Gray"

   399 ; Row header (equally space in y2)

   400   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   401   x3    = x1                         ; same as x1

   402   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   403   text3 = row_header

   404   res3               = True

   405   res3@txFontHeightF = 0.010

   406   res3@gsFillColor   = "Gray"

   408 ; Main table body

   409   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   410   x4    = x2                      ; Start and end x

   411   y4    = y3                      ; Start and end Y

   412   text4 = new((/nrow,ncol/),string)

   414   color_fill4           = new((/nrow,ncol/),string)

   415   color_fill4           = "white"

   416 ; color_fill4(:,ncol-1) = "grey"

   417 ; color_fill4(nrow-1,:) = "green"

   419   res4               = True       ; Set up resource list

   420 ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

   421   res4@txFontHeightF = 0.015

   422   res4@gsFillColor   = color_fill4

   424   delete (color_fill4)

   425 ;-------------------------------------------------------------------

   426 ; for table value

   428   do n = 0,nrow-1

   429      text4(n,0) = sprintf("%5.2f", y81(n+41))

   430      text4(n,1) = sprintf("%5.2f", x81(n+41))

   431      text4(n,2) = sprintf("%5.2f", npp81(n+41))

   432      text4(n,3) = sprintf("%5.2f", rain81(n+41))

   433   end do

   434 ;-------------------------------------------------------------------

   436   plot_name = "table_site81_ob2"

   438   wks = gsn_open_wks (plot_type,plot_name)

   439 ;------------------------------------------

   440 ; for table title

   442   gRes               = True

   443   gRes@txFontHeightF = 0.02

   444 ; gRes@txAngleF      = 90

   446   title_text = "Observation at 81 Sites (2)"

   448   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   449 ;------------------------------------------

   451   gsn_table(wks,ncr1,x1,y1,text1,res1)

   452   gsn_table(wks,ncr2,x2,y2,text2,res2)

   453   gsn_table(wks,ncr3,x3,y3,text3,res3)

   454   gsn_table(wks,ncr4,x4,y4,text4,res4)

   456   frame(wks)

   457   clear (wks)

   459   delete (row_header)

   460   delete (res3)

   461   delete (ncr3)

   462   delete (text3)

   463   delete (res4)

   464   delete (ncr4)

   465   delete (text4)

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

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

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

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

   472 ;**************************************************************************

   473 ;(A)-3 scatter plot, ob 933

   474 ;**************************************************************************

   476  plot_name = "scatter_ob_933"

   477  title     = "Observed 933 sites"

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

   480  res                   = True                  ; plot mods desired

   481  res@tiMainString      = title                 ; add title

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

   483  res@xyMarkers         =  16                   ; choose type of marker

   484  res@xyMarkerColor     = "red"                 ; Marker color

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

   486  res@tmLabelAutoStride = True                  ; nice tick mark labels

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

   489  frame(wks)

   490  clear (wks)

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

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

   494  system("rm "+plot_name+"-*."+plot_type_new)

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

   497 ;**************************************************************************

   498 ;(A)-4 scatter plot, model 933

   499 ;**************************************************************************

   501  plot_name = "scatter_model_933"

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

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

   505  res                   = True                  ; plot mods desired

   506  res@tiMainString      = title                 ; add title

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

   508  res@xyMarkers         =  16                   ; choose type of marker

   509  res@xyMarkerColor     = "red"                 ; Marker color

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

   511  res@tmLabelAutoStride = True                  ; nice tick mark labels

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

   514  frame(wks)

   515  clear (wks)

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

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

   519  system("rm "+plot_name+"-*."+plot_type_new)

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

   522 ;***************************************************************************

   523 ;(A)-5 scatter plot, model vs ob 81

   524 ;***************************************************************************

   526  plot_name = "scatter_model_vs_ob_81"

   527  title     = model_name +" vs ob 81 sites"

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

   530  res                   = True                  ; plot mods desired

   531  res@tiMainString      = title                 ; add title

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

   533  res@xyMarkers         =  16                   ; choose type of marker

   534  res@xyMarkerColor     = "red"                 ; Marker color

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

   536  res@tmLabelAutoStride = True                  ; nice tick mark labels

   538  res@gsnDraw           = False

   539  res@gsnFrame          = False                 ; don't advance frame yet

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

   541 ;compute correlation coef. and M

   542  ccr81 = esccr(nppmod81,npp81,0)

   543 ;print (ccr81)

   545 ;new eq

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

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

   548  print (M81)

   550  tRes  = True

   551  tRes@txFontHeightF = 0.025

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

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

   556 ;--------------------------------

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

   558 ;-------------------------------

   559 ; add polyline

   561  dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)

   562 ;-------------------------------

   563  draw (plot)

   564  frame(wks)

   565  clear (wks)

   566  delete (dum)

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

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

   570  system("rm "+plot_name+"-*."+plot_type_new)

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

   573 ;***************************************************************************

   574 ;(A)-6 scatter plot, model vs ob 933

   575 ;***************************************************************************

   577  plot_name = "scatter_model_vs_ob_933"

   578  title     = model_name +" vs ob 933 sites"

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

   581  res                   = True                  ; plot mods desired

   582  res@tiMainString      = title                 ; add title

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

   584  res@xyMarkers         =  16                   ; choose type of marker

   585  res@xyMarkerColor     = "red"                 ; Marker color

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

   587  res@tmLabelAutoStride = True                  ; nice tick mark labels

   589  res@gsnDraw           = False

   590  res@gsnFrame          = False                 ; don't advance frame yet

   591 ;-------------------------------

   592 ;compute correlation coef. and M

   593  ccr933 = esccr(nppmod933,npp933,0)

   594 ;print (ccr933)

   596 ;new eq

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

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

   599  print (M933)

   601  tRes  = True

   602  tRes@txFontHeightF = 0.025

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

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

   607 ;--------------------------------

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

   609 ;-------------------------------

   610 ; add polyline

   612  dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)

   613 ;-------------------------------

   614  draw (plot)

   615  frame(wks)

   616  clear (wks)

   617  delete (dum)

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

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

   621  system("rm "+plot_name+"-*."+plot_type_new)

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

   624 ;*******************************************************************

   625 ;(A)-7 table of site81 (1) (the first 41 sites), model vs ob

   626 ;*******************************************************************

   628   table_length = 0.95

   630 ; table header name

   631   table_header_name = "Site ID"

   633 ; row (not including header row)

   634   nrow       = 41

   635   row_header = new ((/nrow/),string )

   636   row_header(0:nrow-1) = id81(0:nrow-1)

   638 ; Table header

   639   ncr1  = (/1,1/)               ; 1 row, 1 column

   640   x1    = (/0.005,0.15/)        ; Start and end X

   641   y1    = (/0.85 ,0.95/)        ; Start and end Y

   642   text1 = table_header_name

   643   res1               = True

   644   res1@txFontHeightF = 0.015

   645   res1@gsFillColor   = "CornFlowerBlue"

   647 ; Column header1 (equally space in x2)

   649   delete (col_header)

   650   delete (ncr2)

   651   delete (text2)

   653   col_header = (/"Latitude","Longitude"/)

   654   ncol = dimsizes(col_header)

   656   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   657   x2    = (/x1(1),0.35/)        ; start from end of x1

   658   y2    = y1                    ; same as y1

   659   text2 = col_header

   660   res2               = True

   661   res2@txFontHeightF = 0.015

   662   res2@gsFillColor   = "Gray"

   664 ; Column header2 (equally space in x2)

   666   col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)

   667   col_header2 = (/"ob",model_name \

   668                  ,"ob",model_name \

   669                  /)

   671   ncol1 = dimsizes(col_header1)

   672   ncol2 = dimsizes(col_header2)

   674   ncr21  = (/1,ncol1/)            ; 1 rows, 4 columns

   675   x21    = (/x2(1),0.995/)        ; start from end of x2

   676   yhalf  = (y1(0)+y1(1))*0.5

   677   y21    = (/yhalf,y1(1)/)        ; top half of y1

   678   text21 = col_header1

   679   res21               = True

   680   res21@txFontHeightF = 0.015

   681   res21@gsFillColor   = "Gray"

   683   ncr22  = (/1,ncol2/)            ; 1 rows, 12 columns

   684   x22    = x21                    ; start from end of x1

   685   y22    = (/y1(0),yhalf/)        ; bottom half of y1

   686   text22 = col_header2

   687   res22               = True

   688   res22@txFontHeightF = 0.015

   689   res22@gsFillColor   = "Gray"

   691 ; Row header (equally space in y2)

   692   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   693   x3    = x1                         ; same as x1

   694   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   695   text3 = row_header

   696   res3               = True

   697   res3@txFontHeightF = 0.010

   698   res3@gsFillColor   = "Gray"

   700 ; Main table body-1

   701   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   702   x4    = x2                      ; Start and end x

   703   y4    = y3                      ; Start and end Y

   704   text4 = new((/nrow,ncol/),string)

   706   color_fill4           = new((/nrow,ncol/),string)

   707   color_fill4           = "white"

   708 ; color_fill4(:,ncol-1) = "grey"

   709 ; color_fill4(nrow-1,:) = "green"

   711   res4               = True       ; Set up resource list

   712 ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

   713   res4@txFontHeightF = 0.015

   714   res4@gsFillColor   = color_fill4

   716   delete (color_fill4)

   718 ; Main table body-2

   719   ncr42  = (/nrow,ncol2/)           ; nrow rows, ncol columns

   720   x42    = x21                      ; Start and end x

   721   y42    = y3                       ; Start and end Y

   722   text42 = new((/nrow,ncol2/),string)

   724   color_fill42           = new((/nrow,ncol2/),string)

   725   color_fill42           = "white"

   726 ; color_fill42(:,ncol-1) = "grey"

   727 ; color_fill42(nrow-1,:) = "green"

   729   res42               = True       ; Set up resource list

   730 ; res42@gsnDebug      = True       ; Useful to print NDC row,col values used.

   731   res42@txFontHeightF = 0.015

   732   res42@gsFillColor   = color_fill42

   734   delete (color_fill42)

   735 ;-------------------------------------------------------------------

   736 ; for table value

   738   do n = 0,nrow-1

   739      text4(n,0) = sprintf("%5.2f", y81(n))

   740      text4(n,1) = sprintf("%5.2f", x81(n))

   741   end do

   743   do n = 0,nrow-1

   744      text42(n,0) = sprintf("%5.2f", npp81(n))

   745      text42(n,1) = sprintf("%5.2f", nppmod81(n))

   746      text42(n,2) = sprintf("%5.2f", rain81(n))

   747      text42(n,3) = sprintf("%5.2f", rainmod81(n))

   748   end do

   749 ;---------------------------------------------------------------------------

   751   plot_name = "table_site81_model_vs_ob1"

   753   wks = gsn_open_wks (plot_type,plot_name)

   754 ;------------------------------------------

   755 ; for table title

   757   gRes               = True

   758   gRes@txFontHeightF = 0.02

   759 ; gRes@txAngleF      = 90

   761   title_text = "Model vs Observation at 81 Sites (1)"

   763   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   764 ;------------------------------------------

   766   gsn_table(wks,ncr1,x1,y1,text1,res1)

   767   gsn_table(wks,ncr2,x2,y2,text2,res2)

   768   gsn_table(wks,ncr21,x21,y21,text21,res21)

   769   gsn_table(wks,ncr22,x22,y22,text22,res22)

   770   gsn_table(wks,ncr3,x3,y3,text3,res3)

   771   gsn_table(wks,ncr4,x4,y4,text4,res4)

   772   gsn_table(wks,ncr42,x42,y42,text42,res42)

   774   frame(wks)

   775   clear (wks)

   777   delete (col_header)

   778   delete (row_header)

   779   delete (res1)

   780   delete (ncr1)

   781   delete (text1)

   782   delete (res2)

   783   delete (ncr2)

   784   delete (text2)

   785   delete (res21)

   786   delete (ncr21)

   787   delete (text21)

   788   delete (res22)

   789   delete (ncr22)

   790   delete (text22)

   791   delete (res3)

   792   delete (ncr3)

   793   delete (text3)

   794   delete (res4)

   795   delete (ncr4)

   796   delete (text4)

   797   delete (res42)

   798   delete (ncr42)

   799   delete (text42)

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

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

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

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

   806 ;*******************************************************************

   807 ;(A)-8 table of site81 (2) (the last 40 sites), model vs ob

   808 ;*******************************************************************

   810   table_length = 0.95

   812 ; table header name

   813   table_header_name = "Site ID"

   815 ; row (not including header row)

   816   nrow       = 40

   817   row_header = new ((/nrow/),string )

   818   row_header(0:nrow-1) = id81(41:41+nrow-1)

   820 ; Table header

   821   ncr1  = (/1,1/)               ; 1 row, 1 column

   822   x1    = (/0.005,0.15/)        ; Start and end X

   823   y1    = (/0.85 ,0.95/)        ; Start and end Y

   824   text1 = table_header_name

   825   res1               = True

   826   res1@txFontHeightF = 0.015

   827   res1@gsFillColor   = "CornFlowerBlue"

   829 ; Column header1 (equally space in x2)

   831   col_header = (/"Latitude","Longitude"/)

   832   ncol = dimsizes(col_header)

   834   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   835   x2    = (/x1(1),0.35/)        ; start from end of x1

   836   y2    = y1                    ; same as y1

   837   text2 = col_header

   838   res2               = True

   839   res2@txFontHeightF = 0.015

   840   res2@gsFillColor   = "Gray"

   842 ; Column header2 (equally space in x2)

   844 ; col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)

   845   col_header2 = (/"ob",model_name \

   846                  ,"ob",model_name \

   847                  /)

   849   ncol1 = dimsizes(col_header1)

   850   ncol2 = dimsizes(col_header2)

   852   ncr21  = (/1,ncol1/)            ; 1 rows, 4 columns

   853   x21    = (/x2(1),0.995/)        ; start from end of x2

   854   yhalf  = (y1(0)+y1(1))*0.5

   855   y21    = (/yhalf,y1(1)/)        ; top half of y1

   856   text21 = col_header1

   857   res21               = True

   858   res21@txFontHeightF = 0.015

   859   res21@gsFillColor   = "Gray"

   861   ncr22  = (/1,ncol2/)            ; 1 rows, 12 columns

   862   x22    = x21                    ; start from end of x1

   863   y22    = (/y1(0),yhalf/)        ; bottom half of y1

   864   text22 = col_header2

   865   res22               = True

   866   res22@txFontHeightF = 0.015

   867   res22@gsFillColor   = "Gray"

   869 ; Row header (equally space in y2)

   870   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   871   x3    = x1                         ; same as x1

   872   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   873   text3 = row_header

   874   res3               = True

   875   res3@txFontHeightF = 0.010

   876   res3@gsFillColor   = "Gray"

   878 ; Main table body-1

   879   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   880   x4    = x2                      ; Start and end x

   881   y4    = y3                      ; Start and end Y

   882   text4 = new((/nrow,ncol/),string)

   884   color_fill4           = new((/nrow,ncol/),string)

   885   color_fill4           = "white"

   886 ; color_fill4(:,ncol-1) = "grey"

   887 ; color_fill4(nrow-1,:) = "green"

   889   res4               = True       ; Set up resource list

   890 ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

   891   res4@txFontHeightF = 0.015

   892   res4@gsFillColor   = color_fill4

   894   delete (color_fill4)

   896 ; Main table body-2

   897   ncr42  = (/nrow,ncol2/)           ; nrow rows, ncol columns

   898   x42    = x21                      ; Start and end x

   899   y42    = y3                       ; Start and end Y

   900   text42 = new((/nrow,ncol2/),string)

   902   color_fill42           = new((/nrow,ncol2/),string)

   903   color_fill42           = "white"

   904 ; color_fill42(:,ncol-1) = "grey"

   905 ; color_fill42(nrow-1,:) = "green"

   907   res42               = True       ; Set up resource list

   908 ; res42@gsnDebug      = True       ; Useful to print NDC row,col values used.

   909   res42@txFontHeightF = 0.015

   910   res42@gsFillColor   = color_fill42

   912   delete (color_fill42)

   913 ;-------------------------------------------------------------------

   914 ; for table value

   916   do n = 0,nrow-1

   917      text4(n,0) = sprintf("%5.2f", y81(n+41))

   918      text4(n,1) = sprintf("%5.2f", x81(n+41))

   919   end do

   921   do n = 0,nrow-1

   922      text42(n,0) = sprintf("%5.2f", npp81(n+41))

   923      text42(n,1) = sprintf("%5.2f", nppmod81(n+41))

   924      text42(n,2) = sprintf("%5.2f", rain81(n+41))

   925      text42(n,3) = sprintf("%5.2f", rainmod81(n+41))

   926   end do

   927 ;-------------------------------------------------------------------

   929   plot_name = "table_site81_model_vs_ob2"

   931   wks = gsn_open_wks (plot_type,plot_name)

   932 ;------------------------------------------

   933 ; for table title

   935   gRes               = True

   936   gRes@txFontHeightF = 0.02

   937 ; gRes@txAngleF      = 90

   939   title_text = "Model vs Observation at 81 Sites (2)"

   941   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   942 ;------------------------------------------

   944   gsn_table(wks,ncr1,x1,y1,text1,res1)

   945   gsn_table(wks,ncr2,x2,y2,text2,res2)

   946   gsn_table(wks,ncr21,x21,y21,text21,res21)

   947   gsn_table(wks,ncr22,x22,y22,text22,res22)

   948   gsn_table(wks,ncr3,x3,y3,text3,res3)

   949   gsn_table(wks,ncr4,x4,y4,text4,res4)

   950   gsn_table(wks,ncr42,x42,y42,text42,res42)

   952   frame(wks)

   953   clear (wks)

   955   delete (col_header)

   956   delete (row_header)

   957   delete (res1)

   958   delete (ncr1)

   959   delete (text1)

   960   delete (res2)

   961   delete (ncr2)

   962   delete (text2)

   963   delete (res21)

   964   delete (ncr21)

   965   delete (text21)

   966   delete (res22)

   967   delete (ncr22)

   968   delete (text22)

   969   delete (res3)

   970   delete (ncr3)

   971   delete (text3)

   972   delete (res4)

   973   delete (ncr4)

   974   delete (text4)

   975   delete (res42)

   976   delete (ncr42)

   977   delete (text42)

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

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

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

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

   984 ;**************************************************************************

   985 ;(B) histogram 81

   986 ;**************************************************************************

   987 ;get data

   988   RAIN1_1D = ndtooned(rain81)

   989   RAIN2_1D = ndtooned(rainmod81)

   990   NPP1_1D  = ndtooned(npp81)

   991   NPP2_1D  = ndtooned(nppmod81)

   993 ; number of bin

   994   nx = 8

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

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

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

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

  1000   dx4          = new((/1/),float)

  1002   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

  1003          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

  1005 ;----------------------------------------

  1006 ;compute correlation coeff and M score

  1008  u = yvalues(0,:)

  1009  v = yvalues(1,:)

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

  1012  uu = u(good)

  1013  vv = v(good)

  1015  ccr81h = esccr(uu,vv,0)

  1016 ;print (ccr81h)

  1018 ;new eq

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

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

  1021  print (M81h)

  1023  delete (u)

  1024  delete (v)

  1025  delete (uu)

  1026  delete (vv)

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

  1028 ; histogram res

  1029   resh                = True

  1030   resh@gsnMaximize    = True

  1031   resh@gsnDraw        = False

  1032   resh@gsnFrame       = False

  1033   resh@xyMarkLineMode = "Markers"

  1034   resh@xyMarkerSizeF  = 0.014

  1035   resh@xyMarker       = 16

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

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

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

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

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

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

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

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

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

  1048   lnres = True

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

  1051 ;**************************************************************************

  1052 ;(B)-1 histogram plot, ob 81 site

  1053 ;**************************************************************************

  1055   plot_name = "histogram_ob_81"

  1056   title     = "Observed 81 site"

  1057   resh@tiMainString  = title

  1059   wks   = gsn_open_wks (plot_type,plot_name)

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

  1063 ;-------------------------------

  1064 ;Attach the vertical bar and the horizontal cap line

  1066   delete (x1)

  1067   delete (x2)

  1068   delete (y1)

  1069   delete (y2)

  1071   do nd=0,0

  1072     lnres@gsLineColor = line_colors(nd)

  1073     do i=0,nx-1

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

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

  1077 ;

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

  1079 ;

  1080         x1 = xvalues(nd,i)

  1081         y1 = yvalues(nd,i)

  1082         y2 = mn_yvalues(nd,i)

  1083         plx = (/x1,x1/)

  1084         ply = (/y1,y2/)

  1085         min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1087         y2 = mx_yvalues(nd,i)

  1088         plx = (/x1,x1/)

  1089         ply = (/y1,y2/)

  1090         max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1091 ;

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

  1093 ;

  1094         x1 = xvalues(nd,i) - dx4

  1095         x2 = xvalues(nd,i) + dx4

  1096         y1 = mn_yvalues(nd,i)

  1097         plx = (/x1,x2/)

  1098         ply = (/y1,y1/)

  1099         min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1101         y1 = mx_yvalues(nd,i)

  1102         plx = (/x1,x2/)

  1103         ply = (/y1,y1/)

  1104         max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1105       end if

  1106     end do

  1107   end do

  1109   draw(xy)

  1110   frame(wks)

  1111   clear (wks)

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

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

  1115  system("rm "+plot_name+"-*."+plot_type_new)

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

  1118 ;****************************************************************************

  1119 ;(B)-2 histogram plot, model vs ob 81 site

  1120 ;****************************************************************************

  1122   plot_name = "histogram_model_vs_ob_81"

  1123   title     = model_name+ " vs Observed 81 site"

  1124   resh@tiMainString  = title

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

  1128 ;-----------------------------

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

  1130 ; vertical lines going through the markers.

  1132   resh@pmLegendDisplayMode    = "Always"

  1133 ; resh@pmLegendWidthF         = 0.1

  1134   resh@pmLegendWidthF         = 0.08

  1135   resh@pmLegendHeightF        = 0.05

  1136   resh@pmLegendOrthogonalPosF = -1.17

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

  1138 ; resh@pmLegendParallelPosF   =  0.18

  1139   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

  1141 ; resh@lgPerimOn              = False

  1142   resh@lgLabelFontHeightF     = 0.015

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

  1144 ;-----------------------------

  1145   tRes  = True

  1146   tRes@txFontHeightF = 0.025

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

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

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

  1153 ;-------------------------------

  1154 ;Attach the vertical bar and the horizontal cap line

  1156   do nd=0,1

  1157     lnres@gsLineColor = line_colors(nd)

  1158     do i=0,nx-1

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

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

  1162 ;

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

  1164 ;

  1165         x1 = xvalues(nd,i)

  1166         y1 = yvalues(nd,i)

  1167         y2 = mn_yvalues(nd,i)

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

  1170         y2 = mx_yvalues(nd,i)

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

  1172 ;

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

  1174 ;

  1175         x1 = xvalues(nd,i) - dx4

  1176         x2 = xvalues(nd,i) + dx4

  1177         y1 = mn_yvalues(nd,i)

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

  1180         y1 = mx_yvalues(nd,i)

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

  1182       end if

  1183     end do

  1184   end do

  1186   draw(xy)

  1187   frame(wks)

  1188   clear(wks)

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

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

  1192  system("rm "+plot_name+"-*."+plot_type_new)

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

  1195  delete (RAIN1_1D)

  1196  delete (RAIN2_1D)

  1197  delete (NPP1_1D)

  1198  delete (NPP2_1D)

  1199 ;delete (range)

  1200  delete (xvalues)

  1201  delete (yvalues)

  1202  delete (mn_yvalues)

  1203  delete (mx_yvalues)

  1204  delete (good)

  1205  delete (max_bar)

  1206  delete (min_bar)

  1207  delete (max_cap)

  1208  delete (min_cap)

  1210 ;**************************************************************************

  1211 ;(B) histogram 933

  1212 ;**************************************************************************

  1214 ; get data

  1215   RAIN1_1D = ndtooned(rain933)

  1216   RAIN2_1D = ndtooned(rainmod933)

  1217   NPP1_1D  = ndtooned(npp933)

  1218   NPP2_1D  = ndtooned(nppmod933)

  1220 ; number of bin

  1221   nx = 10

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

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

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

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

  1227   dx4          = new((/1/),float)

  1229   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

  1230          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

  1232 ;----------------------------------------

  1233 ;compute correlation coeff and M score

  1235  u = yvalues(0,:)

  1236  v = yvalues(1,:)

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

  1239  uu = u(good)

  1240  vv = v(good)

  1242  ccr933h = esccr(uu,vv,0)

  1243 ;print (ccr933h)

  1245 ;new eq

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

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

  1248  print (M933h)

  1250  delete (u)

  1251  delete (v)

  1252  delete (uu)

  1253  delete (vv)

  1254 ;----------------------------------------------------------------------

  1255 ; histogram res

  1256   delete (resh)

  1257   resh                = True

  1258   resh@gsnMaximize    = True

  1259   resh@gsnDraw        = False

  1260   resh@gsnFrame       = False

  1261   resh@xyMarkLineMode = "Markers"

  1262   resh@xyMarkerSizeF  = 0.014

  1263   resh@xyMarker       = 16

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

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

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

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

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

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

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

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

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

  1276   lnres = True

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

  1279 ;**************************************************************************

  1280 ;(B)-3 histogram plot, ob 933 site

  1281 ;**************************************************************************

  1283   plot_name = "histogram_ob_933"

  1284   title     = "Observed 933 site"

  1285   resh@tiMainString  = title

  1287   wks   = gsn_open_wks (plot_type,plot_name)

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

  1291 ;-------------------------------

  1292 ;Attach the vertical bar and the horizontal cap line

  1294   do nd=0,0

  1295     lnres@gsLineColor = line_colors(nd)

  1296     do i=0,nx-1

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

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

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

  1303         x1 = xvalues(nd,i)

  1304         y1 = yvalues(nd,i)

  1305         y2 = mn_yvalues(nd,i)

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

  1308         y2 = mx_yvalues(nd,i)

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

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

  1313         x1 = xvalues(nd,i) - dx4

  1314         x2 = xvalues(nd,i) + dx4

  1315         y1 = mn_yvalues(nd,i)

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

  1318         y1 = mx_yvalues(nd,i)

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

  1320       end if

  1321     end do

  1322   end do

  1324   draw(xy)

  1325   frame(wks)

  1326   delete (xy)

  1327   clear (wks)

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

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

  1331  system("rm "+plot_name+"-*."+plot_type_new)

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

  1334 ;****************************************************************************

  1335 ;(B)-4 histogram plot, model vs ob 933 site

  1336 ;****************************************************************************

  1338   plot_name = "histogram_model_vs_ob_933"

  1339   title     = model_name+ " vs Observed 933 site"

  1340   resh@tiMainString  = title

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

  1344 ;-----------------------------

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

  1346 ; vertical lines going through the markers.

  1348   resh@pmLegendDisplayMode    = "Always"

  1349 ; resh@pmLegendWidthF         = 0.1

  1350   resh@pmLegendWidthF         = 0.08

  1351   resh@pmLegendHeightF        = 0.05

  1352   resh@pmLegendOrthogonalPosF = -1.17

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

  1354 ; resh@pmLegendParallelPosF   =  0.18

  1355   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

  1357 ; resh@lgPerimOn              = False

  1358   resh@lgLabelFontHeightF     = 0.015

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

  1360 ;-----------------------------

  1361   tRes  = True

  1362   tRes@txFontHeightF = 0.025

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

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

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

  1369 ;-------------------------------

  1370 ;Attach the vertical bar and the horizontal cap line

  1372   do nd=0,1

  1373     lnres@gsLineColor = line_colors(nd)

  1374     do i=0,nx-1

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

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

  1378 ;

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

  1380 ;

  1381         x1 = xvalues(nd,i)

  1382         y1 = yvalues(nd,i)

  1383         y2 = mn_yvalues(nd,i)

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

  1386         y2 = mx_yvalues(nd,i)

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

  1388 ;

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

  1390 ;

  1391         x1 = xvalues(nd,i) - dx4

  1392         x2 = xvalues(nd,i) + dx4

  1393         y1 = mn_yvalues(nd,i)

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

  1396         y1 = mx_yvalues(nd,i)

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

  1398       end if

  1399     end do

  1400   end do

  1402   draw(xy)

  1403   frame(wks)

  1404   delete(xy)

  1405   clear(wks)

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

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

  1409  system("rm "+plot_name+"-*."+plot_type_new)

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

  1412 ;***************************************************************************

  1413 ;(C) global contour

  1414 ;***************************************************************************

  1416 ;res

  1417   resg                     = True             ; Use plot options

  1418   resg@cnFillOn            = True             ; Turn on color fill

  1419   resg@gsnSpreadColors     = True             ; use full colormap

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

  1421 ; resg@lbLabelAutoStride   = True

  1422   resg@cnLinesOn           = False            ; Turn off contourn lines

  1423   resg@mpFillOn            = False            ; Turn off map fill

  1425   resg@gsnSpreadColors      = True            ; use full colormap

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

  1427   resg@cnMinLevelValF       = 0.              ; Min level

  1428   resg@cnMaxLevelValF       = 2200.           ; Max level

  1429   resg@cnLevelSpacingF      = 200.            ; interval

  1431 ;****************************************************************************

  1432 ;(C)-1 global contour plot, ob

  1433 ;****************************************************************************

  1435   delta = 0.00000000001

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

  1438   plot_name = "global_ob"

  1439   title     = "Observed MODIS MOD 17"

  1440   resg@tiMainString  = title

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

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

  1445   plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)

  1447   frame(wks)

  1448   clear (wks)

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

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

  1452  system("rm "+plot_name+"-*."+plot_type_new)

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

  1455 ;****************************************************************************

  1456 ;(C)-2 global contour plot, model

  1457 ;****************************************************************************

  1459   plot_name = "global_model"

  1460   title     = "Model "+ model_name

  1461   resg@tiMainString  = title

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

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

  1466   plot = gsn_csm_contour_map_ce(wks,nppmod,resg)

  1468   frame(wks)

  1469   clear (wks)

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

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

  1473  system("rm "+plot_name+"-*."+plot_type_new)

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

  1476 ;****************************************************************************

  1477 ;(C)-3 global contour plot, model vs ob

  1478 ;****************************************************************************

  1480   plot_name = "global_model_vs_ob"

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

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

  1485   delete (plot)

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

  1488   resg@gsnFrame             = False          ; Do not draw plot

  1489   resg@gsnDraw              = False          ; Do not advance frame

  1491 ; compute correlation coef and M score

  1493   uu1 = ndtooned(nppmod)

  1494   vv1 = ndtooned(nppglobe)

  1496   delete (good)

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

  1499   ug = uu1(good)

  1500   vg = vv1(good)

  1502   ccrG = esccr(ug,vg,0)

  1503 ; print (ccrG)

  1505   MG   = (ccrG*ccrG)* 5.0

  1506   print (MG)

  1508 ; plot correlation coef

  1510   gRes  = True

  1511   gRes@txFontHeightF = 0.02

  1512   gRes@txAngleF = 90

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

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

  1517 ;--------------------------------------------------------------------

  1518 ;(a) ob

  1520   title     = "Observed MODIS MOD 17"

  1521   resg@tiMainString  = title

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

  1525 ;(b) model

  1527   title     = "Model "+ model_name

  1528   resg@tiMainString  = title

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

  1532 ;(c) model-ob

  1534   zz = nppmod

  1535   zz = nppmod - nppglobe

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

  1538   resg@cnMinLevelValF  = -500           ; Min level

  1539   resg@cnMaxLevelValF  =  500.          ; Max level

  1540   resg@cnLevelSpacingF =  50.           ; interval

  1541   resg@tiMainString    = title

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

  1545   pres                            = True        ; panel plot mods desired

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

  1547                                                 ; indiv. plots in panel

  1548   pres@gsnMaximize                = True        ; fill the page

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

  1552   frame (wks)

  1553   clear (wks)

  1554   delete (plot)

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

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

  1558  system("rm "+plot_name+"-*."+plot_type_new)

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

  1561 ;***************************************************************************

  1562 ;(D)-1 zonal line plot, ob

  1563 ;***************************************************************************

  1565   vv     = zonalAve(nppglobe)

  1566   vv@long_name = nppglobe@long_name

  1568   plot_name = "zonal_ob"

  1569   title     = "Observed MODIS MOD 17"

  1571   resz = True

  1572   resz@tiMainString  = title

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

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

  1578   frame(wks)

  1579   clear (wks)

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

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

  1583  system("rm "+plot_name+"-*."+plot_type_new)

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

  1586 ;****************************************************************************

  1587 ;(D)-2 zonal line plot, model vs ob

  1588 ;****************************************************************************

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

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

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

  1595   s@long_name = nppglobe@long_name

  1596 ;-------------------------------------------

  1597 ; compute correlation coef and M score

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

  1600 ; print (ccrZ)

  1602   MZ   = (ccrZ*ccrZ)* 5.0

  1603   print (MZ)

  1604 ;-------------------------------------------

  1605   plot_name = "zonal_model_vs_ob"

  1606   title     = "Zonal Average"

  1607   resz@tiMainString  = title

  1609   wks = gsn_open_wks (plot_type,plot_name)

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

  1612 ; resz@vpWidthF           = 0.7

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

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

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

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

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

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

  1621   resz@txFontHeightF      = 0.0195            ; change title font heights

  1623 ; Legent

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

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

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

  1627   resz@pmLegendParallelPosF   = .82           ; move units right

  1628   resz@pmLegendOrthogonalPosF = -0.4          ; move units down

  1630   resz@pmLegendWidthF         = 0.10          ; Change width and

  1631   resz@pmLegendHeightF        = 0.10          ; height of legend.

  1632   resz@lgLabelFontHeightF     = .02           ; change font height

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

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

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

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

  1637 ;--------------------------------------------------------------------

  1638   zRes  = True

  1639   zRes@txFontHeightF = 0.025

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

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

  1644 ;--------------------------------------------------------------------

  1646   plot  = gsn_csm_xy (wks,ym,s,resz)

  1648   frame(wks)

  1649   clear (wks)

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

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

  1653  system("rm "+plot_name+"-*."+plot_type_new)

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

  1656 ;***************************************************************************

  1657 ; tar up output plots

  1658 ;***************************************************************************

  1660   temp_name = "temp." + model_name

  1661   system("mkdir -p " + temp_name)

  1662   system("mv *.png " + temp_name)

  1663   system("tar cf "+ temp_name +".tar " + temp_name)

  1665 ;***************************************************************************

  1666 end

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