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

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

     2 ; combine all scatter

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

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

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

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

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

     8 begin

    10  plot_type     = "ps"

    11  plot_type_new = "png"

    13 ;************************************************

    14 ; read in ob data

    15 ;************************************************

    17 ;(A) plot ob scatter_81

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

    20  fil81 = "data.81.nc"

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

    23  x81   = f81->SITE_ID

    24  y81   = f81->TNPP_C

    25  xo81  = f81->LONG_DD

    26  yo81  = f81->LAT_DD

    28  x81@long_name = "SITE_ID"

    29  y81@long_name = "NPP (gC/m2/year)"

    31  plot_name = "scatter_ob_81"

    32  title     = "Observed 81 sites"

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

    35  res                   = True                  ; plot mods desired

    36  res@tiMainString      = title                 ; add title

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

    38  res@xyMarkers         =  16                   ; choose type of marker

    39  res@xyMarkerColor     = "red"                 ; Marker color

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

    41  res@tmLabelAutoStride = True                  ; nice tick mark labels

    43  plot  = gsn_csm_xy (wks,x81,y81,res)              ; create plot

    44  frame(wks)

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

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

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

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

    51  clear (wks)

    53 ;(B) plot ob scatter_933

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

    56  fil933 = "data.933.nc"

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

    59  x933  = f933->SITE_ID

    60  y933  = f933->TNPP_C

    61  xo933 = f933->LONG_DD

    62  yo933 = f933->LAT_DD

    64  x933@long_name = "SITE_ID"

    65  y933@long_name = "NPP (gC/m2/year)"

    67  plot_name = "scatter_ob_933"

    68  title     = "Observed 933 sites"

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

    71  res                   = True                  ; plot mods desired

    72  res@tiMainString      = title                 ; add title

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

    74  res@xyMarkers         =  16                   ; choose type of marker

    75  res@xyMarkerColor     = "red"                 ; Marker color

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

    77  res@tmLabelAutoStride = True                  ; nice tick mark labels

    79  plot  = gsn_csm_xy (wks,x933,y933,res)              ; create plot

    80  frame(wks)

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

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

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

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

    87  clear (wks)

    88 ;************************************************

    89 ; read in model data

    90 ;************************************************

    92 ;(C) plot model scatter_81

    94  model_name = "b30.061n"

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

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

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

   100  ymod = fm->NPP

   101  xm   = fm->lon

   102  ym   = fm->lat

   104  nx = dimsizes(xo81)

   105  do i= 0,nx-1

   106     if (xo81(i) .lt. 0.) then

   107         xo81(i) = xo81(i)+ 360.

   108     end if

   109  end do

   110  print (xo81)

   112  sec_to_year = 86400.*365.

   114  ymod81 = linint2_points(xm,ym,ymod(0,:,:),True,xo81,yo81,0) * sec_to_year

   115  xmod81 = x81

   117  ymod81@long_name = "NPP (gC/m2/year)"

   118  xmod81@long_name = "SITE_ID"

   119 print (ymod81)

   120 print (xmod81)

   122  plot_name = "scatter_model_81"

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

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

   126  res                   = True                  ; plot mods desired

   127  res@tiMainString      = title                 ; add title

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

   129  res@xyMarkers         =  16                   ; choose type of marker

   130  res@xyMarkerColor     = "red"                 ; Marker color

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

   132  res@tmLabelAutoStride = True                  ; nice tick mark labels

   134  plot  = gsn_csm_xy (wks,xmod81,ymod81,res)    ; create plot

   135  frame(wks)

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

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

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

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

   142  clear (wks)

   144 ;(D) plot model scatter_933

   146  model_name = "b30.061n"

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

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

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

   152  ymod = fm->NPP

   153  xm   = fm->lon

   154  ym   = fm->lat

   156  nx = dimsizes(xo933)

   157  do i= 0,nx-1

   158     if (xo933(i) .lt. 0.) then

   159         xo933(i) = xo933(i)+ 360.

   160     end if

   161  end do

   162  print (xo933)

   164  sec_to_year = 86400.*365.

   166  ymod933 = linint2_points(xm,ym,ymod(0,:,:),True,xo933,yo933,0) * sec_to_year

   167  xmod933 = x933

   169  ymod933@long_name = "NPP (gC/m2/year)"

   170  xmod933@long_name = "SITE_ID"

   171 ;print (ymod933)

   172 ;print (xmod933)

   174  plot_name = "scatter_model_933"

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

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

   178  res                   = True                  ; plot mods desired

   179  res@tiMainString      = title                 ; add title

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

   181  res@xyMarkers         =  16                   ; choose type of marker

   182  res@xyMarkerColor     = "red"                 ; Marker color

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

   184  res@tmLabelAutoStride = True                  ; nice tick mark labels

   186  plot  = gsn_csm_xy (wks,xmod933,ymod933,res)    ; create plot

   187  frame(wks)

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

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

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

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

   194  clear (wks)

   196 ;(E) scatter model vs ob 81

   198  ymod81@long_name = "NPP (gC/m2/year)"

   199  y81@long_name    = "NPP (gC/m2/year)"

   200 ;print (ymod81)

   201 ;print (y81)

   203  plot_name = "scatter_model_vs_ob_81"

   204  title     = model_name +" vs ob 81 sites"

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

   207  res                   = True                  ; plot mods desired

   208  res@tiMainString      = title                 ; add title

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

   210  res@xyMarkers         =  16                   ; choose type of marker

   211  res@xyMarkerColor     = "red"                 ; Marker color

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

   213  res@tmLabelAutoStride = True                  ; nice tick mark labels

   214 ;-------------------------------

   215 ;compute correlation coef. and M

   216  ccr81 = esccr(ymod81,y81,0)

   217  print (ccr81)

   219 ;new eq

   220  bias = sum(abs(ymod81-y81)/(ymod81+y81))

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

   222  print (M81)

   224  tRes  = True

   225  tRes@txFontHeightF = 0.025

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

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

   230 ;--------------------------------

   231  plot  = gsn_csm_xy (wks,y81,ymod81,res)       ; create plot

   232  frame(wks)

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

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

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

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

   239  clear (wks)

   241 ;(F) scatter model vs ob 933

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

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

   245 ;print (ymod33)

   246 ;print (y933)

   248  plot_name = "scatter_model_vs_ob_933"

   249  title     = model_name +" vs ob 933 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  ccr933 = esccr(ymod933,y933,0)

   262  print (ccr933)

   264 ;new eq

   265  bias = sum(abs(ymod933-y933)/(ymod933+y933))

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

   267  print (M933)

   269  tRes  = True

   270  tRes@txFontHeightF = 0.025

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

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

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

   276  plot  = gsn_csm_xy (wks,y933,ymod933,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)

   285 end

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