请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 41849|回复: 40

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

  [复制链接]

新浪微博达人勋

发表于 2018-2-12 15:35:12 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
参考
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

wrfout-fixed_height.ncl

3.69 KB, 下载次数: 159, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2021-1-26 17:30:30 | 显示全部楼层
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
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2018-2-12 15:46:04 | 显示全部楼层
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

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-2-12 21:42:32 | 显示全部楼层
先分享,后下载,再学习,谢谢楼主
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-2-12 21:55:56 | 显示全部楼层
先分享,后下载,再学习,谢谢楼主
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-4-12 08:16:03 | 显示全部楼层
谢谢楼主的分享 学习了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-4-19 10:47:47 | 显示全部楼层
请问将经纬度坐标转换成xy坐标的意义在哪里呢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-5-7 22:42:50 | 显示全部楼层
本帖最后由 元宝 于 2019-5-7 23:14 编辑

好棒  非常感谢!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-7-23 21:09:53 | 显示全部楼层
谢谢楼主的分享。为什么我运行的时候,会显示windspdh=spdh(:,x,y)这一行维度不对呢?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-8-21 13:44:18 | 显示全部楼层
不太清楚为什么在使用函数wrf_user_ll_to_ij(a,Longitude,Latitude,res)找到point之后,使用point =point -1呢?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-9-4 07:10:49 | 显示全部楼层
nheight = conform(height,ter,(/0,2,3/)) ; 这句是什么意思呢?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表