NCL学习记录: R RMSE MRSE MAE

7777777lq

于 2023-12-29 16:05:59 发布

阅读量572

点赞数 10

文章标签：学习 linux 笔记

本文链接：https://blog.csdn.net/m0_73705730/article/details/135293286

版权

自己写的检验脚本，是先把wrf转成csv了，因为站点数据缺少15:00的数据，所以尝试用csv对准
（其实可以不用转csv就行！！）

脚本过长，贴上关键代码:

.......

do n = 0, nlines-1
  hours_wrf = where(hours_wrf.eq.15,hours_wrf@_FillValue,\
                    hours_wrf)
end do
  hour_bl = ind(.not.ismissing(hours_wrf))
  nihour_bl = num(.not.ismissing(hours_wrf)) 

if (.not.all(ismissing(hour_bl))) 
    hours = hours_wrf(hour_bl)
    times = TIME(hour_bl)
    lat = LAT(hour_bl)
    lon = LON(hour_bl)
    wspd = stringtofloat(WSPD(hour_bl))
    wdir = stringtofloat(WDIR(hour_bl)) 
  else
    print ("Observations not found: all data are missing")
end if
do i = 0,nlines_ob-1
  month(i)  = where(month(i).eq.month_wrf(0),month_wrf(0),\
                    month@_FillValue) 
end do

if (.not.all(ismissing(month)))
  submon = ind(.not.ismissing(month))
  nimonth = num(.not.ismissing(month))
  sub1 = submon(1)
  sub2 = sub1+nimonth-1   ;submon(nimonth-1)
  ; print(lines_ob(sub2))
else
  print ("Observations not found:"+ month_wrf(0)+"month are all missing")
end if   

  date1 = new(nimonth,integer)
  date1@_FillValue = -999
  date1 = date(sub1:sub2)

do j = 0, nimonth-1
  date1(j)   = where((date1(j).ge.day_wrf(0)).and.(date1(j).lt.day_wrf_end),\
                    date1(j),date1@_FillValue)
end do

  subday = ind(.not.ismissing(date1))
  niday = num(.not.ismissing(date1))

  sub3 = sub1+subday(1)
  print(lines_ob(sub3))
  sub4 = sub3+niday;after day

  delete(month)
  delete(date)
  lines_ob1 = asciiread(dir+filename_ob,-1,"string")
  wd  = stringtofloat(str_get_field(lines_ob1(sub3:sub4),8,delm))  
  wsr = stringtofloat(str_get_field(lines_ob1(sub3:sub4),9,delm))
  month     = stringtoint(str_get_field(lines_ob1(sub3:sub4),2,delm))
  date     = stringtoint(str_get_field(lines_ob1(sub3:sub4),3,delm))
  nwd = dimsizes(wd)

  r = new(dimension_sizes, vartype, parameter)
  rwd = esccr(wd,wspd,0) 
  rwsr = esccr(wsr,wdir,0)
   
  do i = 0,nwd-1
  do while (wd(i).ne.0)
    std_wd(i) = abs((wspd(i)-wd(i))
    std_wsr(i) = abs(wdir(i)-wsr(i))
    stc_wd(i) = std_wd(i)/wd(i)
    stc_wsr(i) = std_wsr(i)/wsr(i)

  end do
end do
    MRSE_wd = sum(stc_wd)
    MRSE_wsr = sum(stc_wsr)

    MAE_wd  = sum(std_wd)
    MAE_wsr  = sum(std_wsr)

    RMSE_wd = sqrt(sum((std_wd(i)*std_wd(i)))
    RMSE_wsr = sqrt(sum((std_wsr(i)*std_wsr(i)))