- 积分
- 1038
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2022-4-5
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 Ronin. 于 2022-10-21 10:21 编辑
我想提取wrfout的固定站点污染物数据与观测值对比,比如说pm10,下面是我观看论坛和官网写的例子,但是一直报错维度不匹配。pm10是四维的,bottom_top是29层,按理说应该是只取第一层,将它变成三维的数据。看帖子应该是使用垂直插值,但是不知道第一层是多少hpa,不知道怎么修改了。
; -------------- 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 in netCDF file and make a loop for all time steps
;********************************************************
in = addfile("/home/gaojian/wrfchem4/WRF/test/em_real/wrfout_d02_2015-12-19_00:00:00","r")
times = wrf_user_list_times(in)
ntimes = dimsizes(times)
lon_sta =111.6507 ;站点经纬度
lat_sta =40.7517
pm10 = new(ntimes,float)
do it = 0, ntimes-1 ;Loop for the time
t = wrf_user_getvar(in, "PM10",it)
;;;返回站点经纬度对应的索引x, y
wres = True
wres@returnInt = True
point = wrf_user_ll_to_ij(in, lon_sta, lat_sta, wres)
point = point-1
x = point(0)
y = point(1)
;;;根据索引读取站点的数据
t_sta = t(y,x)
pm10(it) = t_sta
; T2=t2-273.15
end do
npts=ntimes
fName="pm10.txt"
data = new( npts, "string") ;OUTPUTDATE
print(" Time pm10")
print (times +sprintf("%12.1f", pm10))
data = times +sprintf("%12.1f", pm10)
asciiwrite (fName,data)
end
|
|