- 积分
- 518
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-10-4
- 最后登录
- 1970-1-1
![[哈库拉玛塔塔0000] 粉丝数:112 微博数:10 新浪微博达人勋](source/plugin/sina_login/img/light.png)
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
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()
|
-
|