最原始观测值与模拟值对比计算与绘图

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
file_obs="isop_obs7.csv"
file_sim="oldisopsim.csv"
f_h = asciiread("sitelist.txt",-1,"string")

lines_obs=asciiread(file_obs,-1,"string")
lines_obs@_Fillvalue="-999"
lines_sim=asciiread(file_sim,-1,"string")
delim="," 
R=new(10,float)
MAE=new(10,float)
RMSE=new(10,float)
do i=1,10
time=tofloat(str_get_field(lines_obs,1,delim))
tmp1=str_get_field(lines_obs,i+1,delim)
tmp2=str_get_field(lines_sim,i+1,delim)
allvalue=tofloat(str_get_field(lines_sim,i+1,delim))  ;把所有值都加上
tmp1@_Fillvalue="-999"
do j=0,dimsizes(tmp2)-1
   if (tmp1(j) .eq. "-999") then   
      tmp2(j)=""
      tmp1(j)=""
   end if
end do
var_sim=tofloat(tmp2)  
var_obs=tofloat(tmp1)
;print(var_sim)
;print(var_obs)
iz=ind(.not.ismissing(var_obs))
if (.not.all(ismissing(iz))) then
   R(i-1)=escorc(var_sim(iz),var_obs(iz))
   MAE(i-1)=avg(abs(var_sim(iz)-var_obs(iz)))
   RMSE(i-1)=sqrt(avg((var_sim(iz)-var_obs(iz))^2))
end if
delete (iz)
;print(var_sim)
;print(var_obs)
;print(time)
;print(R)
;print(MAE)
;print(RMSE)
data_all = new((/2,744/),"float")
data_all(0,:)=var_obs
data_all(1,:)=allvalue ;var_sim
wks = gsn_open_wks("x11","gsn_xy"+f_h(i-1))         ; send graphics to PNG file
 
  res               = True  
  res@tiMainString  = f_h(i-1)
  res@tiYAxisString = "Concentration (ppb)"           ; y axis title
  res@tiXAxisString = "Time"                 ; x axis title

  res@xyLineColors      = (/"black","red"/)  ; line colors
  res@xyLineThicknesses = (/4.0,4.0/)        ; line thicknesses
  res@xyDashPatterns    = (/0.0,0.0/)        ; line patterns
  res@xyMarkLineModes   = (/"Lines","Lines"/) ; markers?
  res@xyMarkerColors    = (/"black",    "red"/)     ; color
  res@xyMarkers         = (/16,         2     /)          ; style
  res@tmXBMode          = "Explicit"
  res@tmYLMode          = "Explicit"
  ;res@trYMaxF   = 35
  res@trXMaxF   =744
 ; res@tmYLValues=(/0,20,40,60,80,100,120,140/)
  res@tmXBValues=(/0,96,192,288,384,480,576,672,768/)
  res@tmXBLabels        = (/"07/01","07/05","07/09","07/13","07/17","07/21","07/25","07/29","08/02"/)
 ; res@tmYLLabels        =(/"0","20","40","60","80","100","120","140"/)
  res@tmYLLabelFontHeightF= 0.02
  res@tmXBLabelFontHeightF= 0.02
  res@tiYAxisFontHeightF=0.03
  res@tiXAxisFontHeightF=0.03
  res@vpWidthF = 0.6
  res@vpHeightF= 0.35
  res@pmLegendDisplayMode    = "Always"            ; turn on legend
  res@pmLegendSide           = "Top"               ; Change location of 
  res@pmLegendParallelPosF   = 0.15                  ; move units right
  res@pmLegendOrthogonalPosF = -0.45 ;-0.33                ; move units down
  res@pmLegendWidthF         = 0.15                ; Change width and
  res@pmLegendHeightF        = 0.1                ; height of legend.
  res@lgPerimOn              = False               ; turn off box around
  res@lgLabelFontHeightF     = .03                 ; label font height
  res@xyExplicitLegendLabels = (/"Obs","Sim"/)         ; create explicit labels
 
  plot = gsn_xy(wks,time,data_all,res)          ; Draw an XY plot with 1 curve. 
  txres               = True                     ; text mods desired
  txres@txFontHeightF = 0.03                     ; font smaller. default big
  txres@PerimOn = True
  txres@txBackgroundFillColor = "red"

end do
write_table("stataticsold.txt","w",[/R,MAE,RMSE/],"%8.2f  %8.2f  %8.2f")
end