提取wrfout_10米高度（气压、露点、温度、风向、风速），大佬们帮忙检查下脚本对不对

半岛铁盆

参考自提取wrfout固定高度站点数据（温度、气压、露点、风向、风速））-编程作图-气象家园_气象人自己的家园 (06climate.com)

原代码提取10米风速的时候，输出数据会出现问题（具体见http://bbs.06climate.com/forum.p ... &fromuid=128184）；我在原来的基础上改了代码，能正常输出了，但不知道对不对，大佬们帮忙看下有没有问题。（新人第一次发帖，有什么不对的地方望大佬指正）

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;以下是我修改后的脚本
;提取10米高度（气压、露点、温度、风向、风速）.ncl
; --------------  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-------------------

  dir   = "/home/jina/Build_WRF/WRF输出文档/（2021-09-29）2020-01-01～7号_3级嵌套-方案2/"
  files = systemfunc("ls " + dir + "wrfout_d03_2019*")
  a     = addfiles(files,"r")
  ;a = addfiles("/home/jina/Build_WRF/WRF输出文档/（2021-09-15）2006-04-01～31号/wrfout_d03_2006-04-01_00:00:00" ,"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_ph  = wrf_user_getvar(a,"U10",-1)        ; u averaged to mass points
v_ph  = wrf_user_getvar(a,"V10",-1)        ; v averaged to mass points
tch=wrf_user_getvar(a, "tc",-1)
tdh=wrf_user_getvar(a, "td",-1); td(ntim,nlvl,nlat,nlon)

ph=wrf_user_getvar(a, "pressure",-1)
pbh=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)

;---------------- Select set station lon,lat,fixed height(from surface(m))[CHANGE HERE]-----------------------------
Latitude = 32.067
Longitude = 121.6

;------------------------------------------------------------------
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)
point =point -1
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
nheight = conform(height,ter,(/0,2,3/)) ; assuming height is a 3d array and ter is a 2d array
height = height - nheight

;-------------------------------------------------------------
spdh=sqrt(u_ph^2+v_ph^2)
windspdh=spdh(:,x,y)
dirh=wind_direction(u_ph,v_ph,0)
winddirh=dirh(:,x,y)
tcc=tch(:,0,x,y)
tdc=tdh(:,0,x,y)
phc=ph(:,0,x,y)+pbh(:,0,x,y)  ;units:pa

;  ------ WRITE IN FILE ---------------                    *
npts=ntimes
fName = "结果_提取10米高度（气压、露点、温度、风向、风速）.txt"   ;FILENAME
data  = new( npts, "string")  ;OUTPUTDATE
print("       时间                气压      露点温度      温度         风向        风速 ")
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

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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;以下是修改前的脚本
;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


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

半岛铁盆

自己顶一下

老哥

半岛铁盆 发表于 2021-10-25 10:35
自己顶一下

我用python，可否请教下如何获取对应层的温度，气压，湿度计算公式

老哥

半岛铁盆 发表于 2021-10-25 10:35
自己顶一下

我微信wxid_6jjyxu56pkon22

半岛铁盆

老哥 发表于 2021-11-7 15:37
我用python，可否请教下如何获取对应层的温度，气压，湿度计算公式

我不清楚python怎么求。ncl的话主要就是这两个函数吧wrf_user_getvar、wrf_user_intrp3d，具体怎么计算的你可以看看WRFUserARW.ncl里边的代码。

阿羊出没请注意

你好，我的wrfout数据，没有"pressure","td","cape"等变量，请问WRF输出的变量数是需要设置吗？为什么我没有输出这些变量？

你是最美七月天

请问楼主解决问题了没？

半岛铁盆

阿羊出没请注意 发表于 2022-2-26 20:33
你好，我的wrfout数据，没有"pressure","td","cape"等变量，请问WRF输出的变量数是需要设置吗？为什么我没 ...

我并没有进行特别的设置，输出的wrfout文件中原本就包含气压温度等变量

半岛铁盆

你是最美七月天 发表于 2022-3-20 22:06
请问楼主解决问题了没？

原贴的代码问题我没有解决，但我修改后仅提取10米（不进行向上或向下插值）变量的代码能跑通，没有出现错误，对不对就是另一回事了

木刻

请问我得到的结果维度报错了，该如何解决