将某站点excel里面的观测数据与wrf某格点的数据做plot xy图进行时间上的对比

黄小仙儿

将某站点excel里面的观测数据与wrf某格点的数据做plot xy图进行时间上的对比，横轴为时间
我看ncl 的例子都是直接提取有时间维度的数据直接画图，比如
  u = a->U(0,:,:)，但是现在观测数据是每小时一个数据，没法这么直接提取。
wrf数据已经提取出来了
t2= wrf_user_getvar(a,"T2",-1)
      t2= t2-273.15
      t2@units = "c"
tt2=t2(:,ij(0), ij(1))


但是不知道怎么用观测与模拟做对比，时间变量弄成一致的。

用gsn_csm_xy2画图，求大神指点
还有
plot = gsn_csm_xy2(wks,times,t,p,resL,resR)
中的times,我想以观测数据为准，怎么提取wrf输出的一段时间作为times呢

黄小仙儿

pangzi_xu 发表于 2015-8-18 13:16
我的解决方案是自己设置x轴的时间，只需要wrf的数据和观测数据匹配就好

请问有可以参考的脚本吗？只要一段设置时间的就好。完全没有头绪，谢谢了

黄小仙儿

pangzi_xu 发表于 2015-8-18 16:39
res1D@tmXBMode          = "Explicit"              ; explicit label 这是关键，用好这个就能自己随意 ...

那我想提取一段时间，ntime是72个时次，我只想提取从12开的24个时次呢，ntime怎么提取？

一人。

做XY曲线用NCL貌似不是很好做。可以把模式转出来的数据提取成ascii文件，然后处理起来就比较方便了。

黄小仙儿

pangzi_xu 发表于 2015-8-18 17:33
times_in_file = a->Times
ltime = new(24,string)
  do i=0,23

  ltime(i) = chartostring(times_in_file(i+12,8:12))
这句没怎么看懂。

还有要是我想把几个文件的连起来，比如file1是1号到3号共72个时次，file2是2号到4号共72个时次，file3是3号到5号共72个时次1号24个时次，
我想提取每个file中间的24个时次并且连起来画变量的时序变化图（比如file1提取2号的24个时次，file2提取3号的24个时次。。。。。。）
想了好久不知道怎么实现，大神有什么好方法吗？

黄小仙儿

一人。 发表于 2015-8-18 17:38
做XY曲线用NCL貌似不是很好做。可以把模式转出来的数据提取成ascii文件，然后处理起来就比较方便了。

因为要批量处理，数据都在大型机上，所以不是很方便

黄小仙儿

pangzi_xu 发表于 2015-8-18 21:03
这句是把wrf中的时间变量的第12~36个赋值给ltime数组。如果你把wrf中的时间变量print出来，你会发现8:12 ...

我用自己的方法提取出来了我想要的时间序列，是从wrfout里面提取的。但是最后画图时提示fatal:Argument type mismatch on argument (1) of (gsn_csm_xy) can not coerce
我把时间print出来是这样的，是格式不符合吗？
(701)   2015-05-30_17:00:00
(702)   2015-05-30_18:00:00
(703)   2015-05-30_19:00:00
(704)   2015-05-30_20:00:00
(705)   2015-05-30_21:00:00
(706)   2015-05-30_22:00:00
(707)   2015-05-30_23:00:00
(708)   2015-05-31_00:00:00
(709)   2015-05-31_01:00:00
(710)   2015-05-31_02:00:00
(711)   2015-05-31_03:00:00
(712)   2015-05-31_04:00:00
(713)   2015-05-31_05:00:00
(714)   2015-05-31_06:00:00
(715)   2015-05-31_07:00:00
(716)   2015-05-31_08:00:00
(717)   2015-05-31_09:00:00
(718)   2015-05-31_10:00:00
(719)   2015-05-31_11:00:00
(720)   2015-05-31_12:00:00
(721)   2015-05-31_13:00:00
(722)   2015-05-31_14:00:00
(723)   2015-05-31_15:00:00

黄小仙儿

pangzi_xu 发表于 2015-8-20 15:16
能具体一点吗？你的代码可以发给我看看

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("/public/home/huanglei/wrfdata/wrfout_d03_2015-05-01_12:00:00.nc","r")     

    ij = wrf_user_ll_to_ij(a, 108.58,34.26,True)

    print("lon location is: " +ij(0))

    print("lat location is: " + ij(1))
       
        ; First get the variables we will need
t3=new((/31,24/),float,"No_FillValue")
rain2=new((/31,24/),float,"No_FillValue")
u_2=new((/31,24/),float,"No_FillValue")
v_2=new((/31,24/),float,"No_FillValue")
u10_2=new((/31,24/),float,"No_FillValue")
v10_2=new((/31,24/),float,"No_FillValue")
rh2=new((/31,24/),float,"No_FillValue")
times2=new((/31,24/),string,"No_FillValue")
;rain_tot=new(24,float,"No_FillValue")


filename1=systemfunc("ls /public/home/huanglei/wrfdata/wrfout_d03_2015-05*.nc")
         fin1=addfiles(filename1,"r")
         do i=0,30
            
        t2=fin1->T2(0:23,ij(0), ij(1))
                 rain_exp=fin1->RAINNC(0:23,ij(0), ij(1))
                 rain_con=fin1->RAINC(0:23,ij(0), ij(1))
                 rh0= wrf_user_getvar(fin1,"rh",-1)      
                         u_1=fin1->U(0:23,0,ij(0), ij(1))
                         v_1=fin1->V(0:23,0,ij(0), ij(1))
                 u10_1=fin1->U10(0:23,ij(0), ij(1))
                 v10_1=fin1->V10(0:23,ij(0), ij(1))
                 times = wrf_user_getvar(fin1,"times",-1)



                 tt2= t2-273.15
                 tt2@units = "c"
                 rain_tot= rain_exp + rain_con
                         rh1=rh0(0:23,0,ij(0), ij(1))
               times1=times(0:23)
             ; printVarSummary(times1)
        
              ;print(fin1[1]->T2(0:23,ij(0), ij(1)))
             ;  exit
             t3 (i,:)= tt2         
             rain2(i,:)= rain_tot        
             u_2 (i,:)= u_1        
             v_2(i,:)= v_1        
             u10_2 (i,:)= u10_1         
             v10_2 (i,:)= v10_1         
             rh2 (i,:)= rh1
                times2(i,:)= times1
                end do   
; print(time)
              t=ndtooned(t3)
                          rain=ndtooned(rain2)
                          u=ndtooned(u_2)
                          v=ndtooned(v_2)
                          u10=ndtooned(u10_2)
                          v10=ndtooned(v10_2)
                          rh=ndtooned(rh2)
                          time=ndtooned(times2)



              t_d=t(0:724)
                          rain_d=rain(0:724)
                          u_d=u(0:724)
                          v_d=v(0:724)
                          u10_d=u10(0:724)
                          v10_d=v10(0:724)
                          rh_d=rh(0:724)
                          time_d=time(0:724)


     ;  print(t(0:47))
               
           ;printVarSummary(rh)
           ;printVarSummary(u)
           ;printVarSummary(v)
          ; printVarSummary(u10)
           ;printVarSummary(v10)
           ;printVarSummary(time)
       ;printVarSummary(rain)
   ;     print(time)
      ; print(time)
        ;   exit
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;          
; We generate plots, but what kind do we prefer?
  type = "pdf"
  wks = gsn_open_wks(type,"plt_temp_xy")


; Set some basic resources
  res = True

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


ascii_filename = "/public/home/huanglei/data/ob.txt"
  seismic = asciiread(ascii_filename,(/725,6/),"float")   
  
     time_o= seismic(:,1)  ; Column 1 of file contains X values.
      p_o = seismic(:,2)  ; Column 2 of file contains Y values.
     t2_o= seismic(:,3)  ; Column 3 of file contains Z values.
     ;rh_o= seismic(:,4)  
     wdir_1= seismic(:,4)  
     wspd_1 = seismic(:,5)  
     ;wdir10= seismic(:,7)
    ; wspd10= seismic(:,8)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     t2_oo=t2_o*0.1
         p_oo=p_o *0.1
         wspd=wspd_1*0.1
         rad = 4.0*atan(1.0)/180.
     u_o = -wspd*sin(rad*wdir_1)
     v_o = -wspd*cos(rad*wdir_1)
         
;         u10_o = -wspd10*sin(rad*wdir)
  ;   v10_o = -wspd10*cos(rad*wdir)
  


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
data=new((/2,725/),float)
data(0,:)=t2_oo(:)
data(1,:)=t_d(:)

  print(data)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




res@xyLineThicknesses = (/2.0,2.0/)               ; make 2nd lines thicker
res@gsnYRefLine           = (/-2.0,2.0/)
res@gsnYRefLineThicknesses=(/1.0,1.0/)
;res@xyLineColors      = (/"black","black"/)          ; change line color
res@pmLegendDisplayMode    = "always"              ; turn on legend
res@pmLegendSide           = "Top"                 ; Change location of
res@pmLegendParallelPosF   = .75                 ; move units right
res@pmLegendOrthogonalPosF = -0.25                  ; more neg = down
res@gsnLeftString        ="(a)Temprature"
res@pmLegendWidthF         = 0.10                  ; Change width and
res@pmLegendHeightF        = 0.10                  ; height of legend.
res@lgLabelFontHeightF     = 0.02                   ; change font height
res@lgPerimOn              = False
res@xyExplicitLegendLabels = (/"UB","UF"/)
;res@lgOrientation          ="horizontal"
; res@trXMinF                =1979
; res@trXMaxF                =2011
res@gsnDraw               =False
res@gsnFrame              = False

  
  plot = gsn_csm_xy(wks,time,data,res)

end