提取wrfout固定高度站点数据（温度、气压、露点、风向、风速））

shrubalthea

参考
http://bbs.06climate.com/forum.p ... mp;extra=#pid802903
修改问题后的程序
对齐格式方法：
用atom打开，【另存为】txt
程序内容：
;from wrfout to one station fixed height meteorological data
;history 2018/2/6
;version 1.0
;DEFINE HEIGHT
;variables description reference
;ncl_filedump wrfout_d02_2015-01-01_00_00_00.nc
;http://www.ncl.ucar.edu/Document/Functions/WRF_arw/wrf_user_getvar.shtml

; --------------  LOAD FUNCTIONS AND PROCEDURES ----------------
;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"

begin
; ------------------read wrfout file-------------------
a = addfiles("/home/wrf/WRFV3/run/wrfout_d02_2015-01-01_00_00_00.nc" ,"r")
;a =addfile("/home/RTNWP/wrfout_d02_2014-06-01_00_00_00"+".nc", "r")

; Process all the time steps
times = wrf_user_list_times(a)             ; get all times in the file
ntimes = dimsizes(times)            ; number of times in the file
; -----------read Wind\temperature\dew point temperature\pressure\Height
u  = wrf_user_getvar(a,"ua",-1)        ; u averaged to mass points
v  = wrf_user_getvar(a,"va",-1)        ; v averaged to mass points
tc=wrf_user_getvar(a, "tc",-1)
td=wrf_user_getvar(a, "td",-1); td(ntim,nlvl,nlat,nlon)
;rh=wrf_user_getvar(a, "rh",-1);relative humidity
;psfc=wrf_user_getvar(a, "PSFC",-1) ;surface Pressure
p=wrf_user_getvar(a, "pressure",-1)
pb=wrf_user_getvar(a, "PB",-1)
height  = wrf_user_getvar(a, "z",-1)
ter = wrf_user_getvar(a, "ter",-1) ; ter=HGT(HGT_M, HGT_U, HGT_V)
; printVarSummary (psfc)

;---------------- Select set station lon,lat,fixed height(from surface(m))[CHANGE HERE]-----------------------------
Latitude = 36.05
Longitude = 103.88
hh=100
;------------------------------------------------------------------

res = True
res@returnInt = True                       ; False : return real values, True: return interger values
point = wrf_user_ll_to_ij(a,Longitude,Latitude,res)       ; wrf_user_ll_to_ij(nc_file,lon,lat,opt)
;The function wrf_user_ll_to_ij find the nearest grid point for a specific lat and lon
point =point -1
x = point(0)
y = point(1)
;print("X location is:" + x)
;print("Y location is:" + y)
nheight = conform(height,ter,(/0,2,3/)) ; assuming height is a 3d array and ter is a 2d array
height = height - nheight

;-------------------interpolation----------------------
;wrf_user_intrp3d(三维数组，垂直数组（压力\高度），插值信息，TRUE=从point A点到point B的横截面图;否则为False。)
u_ph  = wrf_user_intrp3d( u,height,"h", hh,0.,False)
v_ph  = wrf_user_intrp3d( v,height,"h", hh,0.,False)
tch   = wrf_user_intrp3d( tc,height,"h", hh,0.,False)
tdh   = wrf_user_intrp3d( td,height,"h", hh,0.,False)
ph   = wrf_user_intrp3d(p,height,"h", hh,0.,False)
pbh  =wrf_user_intrp3d(pb,height,"h", hh,0.,False)
;rhh=wrf_user_intrp3d(rh,height,"h", hh,0.,False)
;-------------------------------------------------------------
spdh= sqrt(u_ph^2 + v_ph^2)
windspdh=spdh(:,x,y)
dirh=wind_direction(u_ph,v_ph,0)
winddirh=dirh(:,x,y)
;print(winddirh)
tcc=tch(:,x,y)
tdc=tdh(:,x,y)
phc=ph(:,x,y)+pbh(:,x,y)  ;units:pa
;psfcp=psfc(:,x,y)
; rhc=rhh(:,x,y)
;  ------ WRITE IN FILE ---------------                    *
npts=ntimes
fName = "heightdata.txt"   ;FILENAME
data  = new( npts, "string")  ;OUTPUTDATE
print("  Time   pressure  dew point temperature  Temperature  Wind_dir_hm      Wind_speed_hm ")
print ( times+sprintf("%12.1f",phc/100)+sprintf("%12.1f",tdc)+ sprintf("%12.1f",tcc)+sprintf("%12.0f",winddirh)+sprintf("%12.1f", windspdh) )
data = times +sprintf("%12.1f",phc/100)+sprintf("%12.1f",tdc)+ sprintf("%12.1f",tcc)+ sprintf("%12.0f",winddirh)+sprintf("%12.1f", windspdh)
asciiwrite (fName,data)
end

caiyjshine

tinalovelife 发表于 2019-8-21 13:44
不太清楚为什么在使用函数wrf_user_ll_to_ij(a,Longitude,Latitude,res)找到point之后，使用point =point - ...

ll to ij是1based，也就是从1开始计数。
ll to xy是0based，也就是从0开始计数。
ncl默认的数组从0开始计数。所以推荐用xy

shrubalthea

CTRL+A,CTRL+C,新建，CTRL+V,另存为xx.txt  即可，不需要用FORTRAN调整格式
格式：                               气压       露点          温度         风向          风速
2015-01-01_12:00:00         801.5       -17.0           0.3           70           3.3
2015-01-01_13:00:00         801.5       -17.8          -0.0           62           4.4
2015-01-01_14:00:00         801.5       -15.5          -1.0           53           4.1

13top

先分享，后下载，再学习，谢谢楼主

13top

先分享，后下载，再学习，谢谢楼主

放逐之刃

谢谢楼主的分享 学习了

chenyf

请问将经纬度坐标转换成xy坐标的意义在哪里呢

元宝

好棒  非常感谢！

羽陌轻寒

谢谢楼主的分享。为什么我运行的时候，会显示windspdh=spdh(:,x,y)这一行维度不对呢？

tinalovelife

不太清楚为什么在使用函数wrf_user_ll_to_ij(a,Longitude,Latitude,res)找到point之后，使用point =point -1呢？

leizhang_zl

nheight = conform(height,ter,(/0,2,3/)) ; 这句是什么意思呢？