爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 3233|回复: 3

wrfout提取的风场数据存为nc格式,不能用MeteoInfoMap打开

[复制链接]

新浪微博达人勋

发表于 2022-3-10 16:21:40 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 哈库拉玛塔塔000 于 2022-3-10 16:23 编辑

编写的以下脚本,生成的nc数据不能用MeteoInfoMap打开,不知是哪里写得不规范,还请大佬指教!
inputPath = "F:\\wrfout_d02_2022-03-10_16_00_00"
f = addfile(inputPath)
time = inputPath.split('_')[2]+inputPath.split('_')[3]+':'+inputPath.split('_')[4]+':'+inputPath.split('_')[5]
#Calculate pressure array
pp = f['P'][:]
psb = f['PB'][:]
ptop = 5000.0
eta = f['ZNU'][0]
latitude = f['XLAT'][:]
longitude = f['XLONG'][:]
nz,ny,nx = pp.shape
nt = 1
eta = eta.reshape(nz,1,1)
pres = (psb * eta + ptop + pp) * 0.01
#Read and de-stagger u and v
u = f['U'][:]
v = f['V'][:]
destag_u = meteo.wrf.destagger(u, -1)
destag_v = meteo.wrf.destagger(v, -2)
y = destag_u.dimvalue(-2)
x = destag_u.dimvalue(-1)
#Interpolate u and v to pressure level
levels = [600, 500, 200]
nz = len(levels)
uwind = meteo.log_interpolate_1d(levels, pres, destag_u, axis=0)
vwind = meteo.log_interpolate_1d(levels, pres, destag_v, axis=0)

#Save netcdf file
outfn = ""
for i in range(0, len(inputPath.split('\\'))):
    if i==len(inputPath.split('\\'))-1:
        outfn = outfn + inputPath.split('\\')
    else:
        outfn = outfn + inputPath.split('\\')+"\\"
outfn = outfn+'.NC'
ncfile = addfile(outfn, 'c')

tidim = ncfile.adddim('Time', nt)
ledim = ncfile.adddim('bottom_top', nz)
sndim = ncfile.adddim('south_north', ny)
wedim = ncfile.adddim('west_east', nx)

ncfile.addgroupattr('ITLE', ' OUTPUT FROM *             PROGRAM:WRF-Chem V4.0.3 MODEL')
ncfile.addgroupattr('START_DATE', time)
ncfile.addgroupattr('SIMULATION_START_DATE', time)
ncfile.addgroupattr('WEST-EAST_GRID_DIMENSION', nx)
ncfile.addgroupattr('SOUTH-NORTH_GRID_DIMENSION', ny)
ncfile.addgroupattr('BOTTOM-TOP_GRID_DIMENSION', nz)
ncfile.addgroupattr('DX', '1000.0f')
ncfile.addgroupattr('DY', '1000.0f')

variables = []
var0 = ncfile.addvar('XLAT', 'float', [tidim, sndim, wedim])    #Latitude
var0.addattr('FieldType', 104)
var0.addattr('description', 'LATITUDE, SOUTH IS NEGATIVE')
var0.addattr('units', 'degree_north')
var0.addattr('coordinates', 'XLONG XLAT')
variables.append(var0)

var1 = ncfile.addvar('XLONG', 'float', [tidim, sndim, wedim])    #Longitude
var1.addattr('FieldType', 104)
var1.addattr('description', 'LONGITUDE, SOUTH IS NEGATIVE')
var1.addattr('units', 'degree_east')
var1.addattr('coordinates', 'XLONG XLAT')
variables.append(var1)

var2 = ncfile.addvar('LEVEL', 'int', [tidim, ledim])    #Level
var2.addattr('FieldType', 104)
var2.addattr('description', 'PRESSURE LEVEL')
var2.addattr('units', 'hPa')
variables.append(var2)

var3 = ncfile.addvar('U', 'float', [tidim, ledim, sndim, wedim])    #U
var3.addattr('FieldType', 104)
var3.addattr('description', 'x-wind component')
var3.addattr('units', 'm s-1')
var3.addattr('coordinates', 'XLONG_U XLAT_U XTIME')
variables.append(var3)

var4 = ncfile.addvar('V', 'float', [tidim, ledim, sndim, wedim])   #V
var4.addattr('FieldType', 104)
var4.addattr('description', 'y-wind component')
var4.addattr('units', 'm s-1')
var4.addattr('coordinates', 'XLONG_V XLAT_V XTIME')
variables.append(var4)

ncfile.create()

ncfile.write(variables[0], latitude.reshape(nt,ny,nx))
ncfile.write(variables[1], longitude.reshape(nt,ny,nx))
ncfile.write(variables[2], array(levels).reshape(nt,nz))
ncfile.write(variables[3], uwind.reshape(nt,nz,ny,nx))
ncfile.write(variables[4], vwind.reshape(nt,nz,ny,nx))

ncfile.flush()
ncfile.close()

错误.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2022-3-10 17:51:14 | 显示全部楼层
用NCO提取的就可以
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-3-15 10:29:29 | 显示全部楼层
已解决,将原数据的全局属性全部赋给生成的数据。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-3-23 14:20:26 | 显示全部楼层
留存,以后需要用
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

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