从wrfout提取指定高度和经纬度的风向和风速到txt文件

physzw

本贴是参考http://bbs.06climate.com/forum.php?mod=viewthread&tid=54024，在脚本里加上了风向的计算公式，可以提取出指定高度和经纬度的风向和风速到txt文件。
我学习NCL时间不长，非常感谢大家的帮助。
    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/csm/contributed.ncl"
    load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

    ; --------------  BEGINING OF NCL SCRIPT ----------------


    begin
    ;********************************************************
    ; read in netCDF file
    ;********************************************************
      
      in   = addfiles("./wrfout_d01_2016-01-01_00:00:00","r")
      ListSetType (in, "cat")  


    ;********************************************************
    ; Process all the time steps
    ;********************************************************



    times = wrf_user_list_times(in)             ; get all times in the file
    ntimes = dimsizes(times)            ; number of times in the file

  
   ; Wind and Height
    u  = wrf_user_getvar(in,"ua",-1)        ; u averaged to mass points
    v  = wrf_user_getvar(in,"va",-1)        ; v averaged to mass points
    height  = wrf_user_getvar(in, "z",-1) ; height is our vertical coordinate
    ter = wrf_user_getvar(in, "ter",-1) ; model terrain height (HGT_M, HGT_U, HGT_V)


    printVarSummary (u)
    printVarSummary (v)
    printVarSummary (height)
    printVarSummary (ter)


    ;************************************************
    ;  - Select lon & lat of the point of interest -
    ;************************************************
    ; - The function wrf_user_ll_to_ij find the nearest grid point for a specific lat and lon

    Latitude = 36.65
    Longitude = 117   

    res = True      
    res@returnInt = True                       ; False : return real values, True: return interger values
    point = wrf_user_ll_to_ij(in,Longitude,Latitude,res)       ; wrf_user_ll_to_ij(nc_file,lon,lat,opt)
     
     x = point(0)
     y = point(1)

     ;print("X location is: " + x)            ; print the value of X at the screen
     ;print("Y location is: " + y)            ; print the value of Y at the screen



    ;*************************************************************************************
    ;  - extract wind, Temperature, Pressure, relative humidity and height coordinates-  *
    ;*************************************************************************************


    ; Conform data to Terrain Height   
       nheight = conform(height,ter,(/0,2,3/)) ; assuming height is a 3d array and ter is a 2d array
       height = height - nheight
   
    ;*******************************************************************************
    ;     - Interpolate wind speed and wind direction at different heights -          *
    ;*******************************************************************************

                    ; Interpolate U,V to 100 Meters
          u_p100  = wrf_user_intrp3d( u,height,"h", 100,0.,False)
          v_p100  = wrf_user_intrp3d( v,height,"h", 100,0.,False)
     
          ; Calculate Wind Speed from Vectors
          spd100 = sqrt(u_p100^2 + v_p100^2)
          printVarSummary(spd100)
          windspd100=spd100(:,x,y)
          dire=180+atan2(u_p100,v_p100)*180/3.14
          wdir100=dire(:,x,y)
         
          printVarSummary(windspd100)


    ;************************************************************
    ;  - Print the variables at the screen -                    *
    ;************************************************************

    npts=ntimes
    fName = "wind100.txt"
    data  = new( npts, "string")
      
     
    print("  Time      Wind_dir_100m      Wind_speed_100m ")
     

        print ( times   + " " +sprintf("%12.2f", wdir100) +"  " +sprintf("%12.2f", windspd100) )                     



       data = times  +" " +sprintf("%12.2f", wdir100) +"  " +sprintf("%12.2f", windspd100)

asciiwrite (fName , data)

end

大气杨娃娃

感谢对前人工作的补充和完善。

羽陌轻寒

感谢楼主分享，为什么我运行的时候会在  windspd100=spd100(:,x,y) 这一行报错呢？

右有枯黄

官网的 wrf_user_ll_to_ij里面有-1的命令
应该是  point = wrf_user_ll_to_ij(in,Longitude,Latitude,res)       ; wrf_user_ll_to_ij(nc_file,lon,lat,opt)
     
     x = point(0)-1
     y = point(1)-1

huasheng

风向风速直接有现成函数

physzw

huasheng 发表于 2019-9-29 10:18
风向风速直接有现成函数

你好，可以详细说明一下吗？谢谢

笑个YY

羽陌轻寒 发表于 2019-7-23 17:25
感谢楼主分享，为什么我运行的时候会在  windspd100=spd100(:,x,y) 这一行报错呢？

请问解决了吗

九十

huasheng 发表于 2019-9-29 10:18
风向风速直接有现成函数

请问为什么我现在使用wind_speed函数会被提示这个函数未被定义呢？请问您有遇到过这种问题吗，我的ncl版本是6.6.2，我看官网上说这个函数6.4.0之后版本都能用，所以应该不是版本的问题，且wind_direction函数也不会报错。

ppprs

九十 发表于 2024-1-24 10:01
请问为什么我现在使用wind_speed函数会被提示这个函数未被定义呢？请问您有遇到过这种问题吗，我的ncl版 ...

您好 请问您是如何提取风速文件的？