# #!/usr/bin/env python
#'''
# Copyright 20152019 ECMWF.
#
# This software is licensed under the terms of the Apache Licence Version 2.0
# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0
#
# In applying this licence, ECMWF does not waive the privileges and immunities
granted to it by
# virtue of its status as an intergovernmental organisation
nor does it submit to any jurisdiction.
#
# **************************************************************************
# Function : compute_geopotentialwind_onspeed_mlheight
#
# Author (date) : Cristian Simarro (0920/1011/2015)
#
#modified: Category : COMPUTATION
#
# OneLineDescXavi Abellan : Computes geopotential on model levels
#
# Description (03/12/2018) - compatibilty with Python 3
Category : COMPUTATION
OneLineDesc : Computes geopotentialthe onu/v model levels.
# components of the wind at certain height
Description : Computes computes the u/v components of the wind Basedat oncertain
code from Nils Wedi, the IFS documentation:
# height. First, it calculates the geopotential of https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation
#each model
level. partOnce III.the Dynamicsrequested andheight numericalis procedures
#found between two model
optimised implementation bylevels, Dominique Lucas.
# the program will vertically interpolate the u/v
ported to Python by Cristiancomponent Simarro
#
# Parametersvalues.
: tq.grib Based on code from Nils Wedi, the IFS documentation:
https://www.ecmwf.int/en/forecasts/documentation-and-support/changes-ecmwf-model/ifs-documentation
grib file with all the levelist of t and q
# part III. Dynamics and numerical procedures
zlnsp.grib optimised implementation by -Dominique gribLucas.
file with levelist 1 for params z and lnsp
# ported to Python by Cristian Simarro
Parameters : -ow outputwind (optional) - name of the output file (default='z_out.grib')
#
# Return Value : output (default='z_out.grib')
# - height in meters you want to know
A fieldset of geopotential on model levels
#
# Dependencies : None
#
# Example Usage :
# python compute_geopotential_on_ml.py tq.grib zlnsp.grib
from numpy import *
import sys,math,os
import argparse
from gribapi import *
def main(t_q,z_lnsp,u_v,out_name,h):
#some checks and information printing
print "Using as input files:\n ",t_q,z_lnsp,u_v
print "The result will be stored in:\n ",out_name
if os.path.exists(out_name):
u/v components
tq.grib - grib file with all the levelist of
os.remove(out_name)
fout = open(out_name,'w')
ftmp = open(t_q)
ftmp.close()
Rdt =and 287.06q
RG = 9.80665
index_keys = ["date","time","shortName","level","step"]
uv.grib values= {}
pv = {}
out = {} - grib file with all the levelists of
gid_out = {}
values_plev = {}
values_lev = {}
h=int(h)
zlnsp = grib_index_new_from_file(z_lnsp,index_keys) u/v
iidtq = grib_index_new_from_file(t_q,index_keys)
iiduv = grib_index_new_from_file(u_v,index_keys)
zlnsp.grib
#we need to get z- andgrib lnspfile fromwith thelevelist first1 levelfor
to do the calculations
counter=0
for date in grib_index_get(zlnsp,'date'):
grib_index_select(zlnsp,'date',date)
grib_index_select(iidtq,'date',date)
grib_index_select(iiduv,'date',date)
params z and lnsp
for time in grib_index_get(zlnsp,'time'):
-o output grib_index_select(zlnsp,'time',time)
(optional) - name of the output file
grib_index_select(iidtq,'time',time)
grib_index_select(iiduv,'time',time)
grib_index_select(zlnsp,'level',1)
grib_index_select(zlnsp,'step',0)
grib_index_select(zlnsp,'shortName','z(default='uv_out_<wind>.grib')
Return Value : output (default='uv_out_<wind>.grib')
gid = grib_new_from_index(zlnsp)
A fieldset the u/v components values for at the specified
#gridType must be gridded, not spectral
height
Dependencies if grib_get(gid,"gridType",str) == "sh": None
Example Usage :
print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral')
sys.exit(1)
compute_wind_speed_height.py tq.grib zlnsp.grib uv.grib -w 100
'''
from __future__ import print_function
import sys
import argparse
import numpy as np
from eccodes import (codes_index_new_from_file, codes_index_get, codes_get,
# surface geopotential
codes_index_select, codes_new_from_index, codes_set,
values["z"] = grib_get_values(gid)
z_h = values["z"]codes_index_add_file, codes_get_array, codes_get_values,
pv = grib_get_array(gid,'pv')
codes_index_release, levelSizeNV = grib_get(gid,'NV',int)/2 -1codes_release, codes_set_values,
grib_release(gid)
codes_write)
R_D = 287.06
R_G = 9.80665
def parse_args():
''' Parse program forarguments step in grib_index_get(iidtq,'step'):using ArgumentParser'''
parser = argparse.ArgumentParser(
z_h = values["z"]
# heigh in geopotential
description='Python tool to calculate the wind at certain height')
parser.add_argument('-w', '--height', required=True, type=int,
my_z = h*RG + z_h
help='height to calculate the wind components')
parser.add_argument('-o', '--output', help='name of the output file')
z_f_prev = z_h parser.add_argument('t_q', metavar='tq.grib', type=str,
help=('grib file with temperature(t) grib_index_select(iidtq,'step',step)and humidity(q)',
grib_index_select(iiduv,'step',step)
'for the for shortName in ["lnsp"]:
model levels'))
parser.add_argument('z_lnsp', metavar='zlnsp.grib', type=str,
grib_index_select(zlnsp,'shortName',shortName)
help=('grib file with geopotential(z) and Logarithm',
grib_index_select(zlnsp,'step',step)
'of surface pressure(lnsp) for the ml=1'))
parser.add_argument('u_v', metavar='uv.grib', type=str,
gid = grib_new_from_index(zlnsp)
if grib_get(gid,"gridType",str) == "sh":
help=('grib file with u and v component of wind for',
print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral')
'the model levels'))
args = sysparser.exitparse_args(1)
if not args.output:
args.output values[shortName] = grib_get_values(gid)= 'uv_out_%sm.grib' % args.height
return args
def main():
'''Main function'''
args pv = gribparse_get_array(gid,'pv'args()
print('Arguments: %s' % ", ".join(
['%s: levelSizeNV = grib_get(gid,'NV',int)/2 -1
grib_release(gid)
%s' % (k, v) for k, v in vars(args).items()]))
fout = open(args.output, 'wb')
index_keys = ['date', 'time', # surface pressure'shortName', 'level', 'step']
idx = codes_index_new_from_file(args.z_lnsp, index_keys)
codes_index_add_file(idx, args.t_q)
if sp = exp(values["lnsp"])
'u_v' in args:
codes_index_add_file(idx, args.u_v)
# iterate date
# get the coefficientsfor fordate computing the pressuresin codes_index_get(idx, 'date'):
codes_index_select(idx, 'date', date)
# how many# levels are we computing?iterate step
for time gribin codes_index_get(idx, 'time'):
codes_index_select(iidtqidx, 'shortNametime',"t" time)
values levelSize=max(grib_index_get(iidtq,"level",int))= get_initial_values(idx, keep_sample=True)
if levelSize != levelSizeNV'height' in args:
values['height'] = args.height
print(sys.argv[0]+' [WARN] total levels should be: '+str(levelSizeNV)+' but it is '+str(levelSize))
values['gh'] = args.height * R_G + values['z']
A = pv[0:levelSize+1]
if 'levelist' in args:
B = pv[levelSize+1:]
Ph_levplusone = A[levelSize] + (B[levelSize]*sp)
values['levelist'] = args.levelist
# Weiterate wantstep toall integratebut upgeopotential intoz thewhich atmosphere,is startingalways atstep the0 ground(an)
for step in # so we start at the lowest level (highest number) and keep
codes_index_get(idx, 'step'):
codes_index_select(idx, 'step', step)
# accumulating the height as we go.
# surface pressure
# See the IFS documentationtry:
# https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation
values['sp'] = get_surface_pressure(idx)
# part III
production_step(idx, values, fout)
# For speed and file I/O, we performexcept theWrongStepError:
computations with numpy vectors instead
if #step of fieldsets.!= '0':
raise
#initializetry:
values for the output
codes_release(values['sample'])
for param in ["u","v"]except KeyError:
pass
codes_index_release(idx)
grib_index_select(iiduv,'level',1)fout.close()
def get_initial_values(idx, keep_sample=False):
'''Get the values of surface z, pv and number of levels '''
gribcodes_index_select(iiduvidx, 'shortNamelevel',param 1)
codes_index_select(idx, 'step', 0)
codes_index_select(idx, 'shortName', 'z')
gid = gribcodes_new_from_index(iiduvidx)
values = {}
# surface geopotential
values['z'] = gid_out[param]=grib_clonecodes_get_values(gid)
values['pv'] = codes_get_array(gid, 'pv')
values['nlevels'] = codes_get(gid, 'NV', int) // 2 grib_release(gid)- 1
check_max_level(idx, values)
if keep_sample:
out[param]=zeros(sp.size)values['sample'] = gid
else:
codes_release(gid)
found = [False for i in range(sp.size)]return values
def check_max_level(idx, values):
'''Make sure we have all the levels required'''
# how formany levlevels in list(reversed(range(1,levelSize+1))):are we computing?
max_level = max(codes_index_get(idx, 'level', int))
if max_level != values['nlevels']:
print('%s [WARN] total levels should #be: select%d thebut levelistit andis retrieve%d' the%
vaules of t and q
(sys.argv[0], # t_level: values for t
values['nlevels'], max_level),
file=sys.stderr)
# q_level: values for q
gribvalues['nlevels'] = max_level
def get_surface_pressure(idx):
'''Get the surface pressure for date-time-step'''
codes_index_select(iidtqidx, 'level',lev 1)
gribcodes_index_select(iidtqidx, 'shortName',"t" 'lnsp')
gid = gribcodes_new_from_index(iidtqidx)
if gid is None:
raise WrongStepError()
t_level = gribif codes_get_values(gid, 'gridType', str) == 'sh':
print('%s [ERROR] fields must be gridded, not spectral' % sys.argv[0],
grib_release(gid)
file=sys.stderr)
sys.exit(1)
grib_index_select(iidtq,'shortName',"q")
# surface pressure
sfc_p = np.exp(codes_get_values(gid))
codes_release(gid)
return gid = grib_new_from_index(iidtq)sfc_p
def get_ph_levs(values, level):
'''Return the presure at a given level and the next'''
qa_levelcoef = grib_get_values(gid)
values['pv'][0:values['nlevels'] + 1]
b_coef = values['pv'][values['nlevels'] + 1:]
ph_lev = a_coef[level - 1] + grib_release(gid)
(b_coef[level - 1] * values['sp'])
ph_levplusone = a_coef[level] + (b_coef[level] * values['sp'])
return ph_lev, # compute moist temperatureph_levplusone
def compute_z_level(idx, lev, values, z_h):
'''Compute z at half & full level for the given level, based on t_level = t_level * (1.+0.609133*q_level)
t/q/sp'''
# select the levelist and retrieve the vaules of t and q
# compute the pressures (on half-levels)t_level: values for t
# q_level: values for q
codes_index_select(idx, 'level', lev)
Ph_lev = A[lev-1] + (B[lev-1] * sp)
codes_index_select(idx, 'shortName', 't')
gid = codes_new_from_index(idx)
t_level = codes_get_values(gid)
if lev == 1:codes_release(gid)
codes_index_select(idx, 'shortName', 'q')
gid = codes_new_from_index(idx)
q_level = codes_get_values(gid)
dlogP = log(Ph_levplusone/0.1)codes_release(gid)
# compute moist temperature
t_level = t_level * (1. + 0.609133 * q_level)
# alphacompute the =pressures log(2on half-levels)
ph_lev, ph_levplusone = get_ph_levs(values, lev)
if lev else== 1:
dlog_p = np.log(ph_levplusone / 0.1)
alpha dlogP = np.log(Ph_levplusone/Ph_lev2)
else:
dlog_p = np.log(ph_levplusone / ph_lev)
alpha dP= 1. - ((ph_lev =/ Ph(ph_levplusone -Ph ph_lev
)) * dlog_p)
t_level = t_level * R_D
# z_f is the geopotential of this full level
alpha # =integrate 1.from -previous ((Ph_lev/dP)*dlogP)
lower) half-level z_h to the
# full level
z_f = z_h TRd = + (t_level *Rd alpha)
# z_h is the geopotential of 'half-levels'
# integrate z_h to next half #level
z_f is the geopotentialz_h of= thisz_h full+ (t_level
* dlog_p)
return z_h, z_f
def production_step(idx, values, fout):
'''Produce u/v interpolated from ML #for integratea from previous (lower) half-level z_h to the full level
given height'''
# We want to integrate up into the atmosphere, starting at the
# ground so we start at the lowest level (highest number) and
z_f = z_h + (TRd*alpha)
# keep accumulating the height as we go.
# See the IFS documentation, part III
# z_hFor isspeed theand geopotential of 'half-levels'
file I/O, we perform the computations with
# numpy vectors instead of fieldsets.
out = {}
# integrate z_h toparams next half level= ('u', 'v')
# we need to get z and lnsp from the first level to do the
z_h=z_h+(TRd*dlogP)
# calculations
z_h = values['z']
# height in geopotential
Ph_levplusonez_f_prev = Phz_levh
# initialize values for the output
for param in params:
# store the result (z_f) in a field and add to the output fieldset
out[param] = np.zeros(values['sp'].size)
found = [False for i in range(values['sp'].size)]
for lev in list(reversed(list(range(1, values['nlevels'] + 1)))):
# (add it to the front, not the end, because we are going 'backwards' z_h, z_f = compute_z_level(idx, lev, values, z_h)
# retrieve u/v params for the # through the fields)current level
for param in params:
for param in ["u","v"]:
codes_index_select(idx, 'level', lev) # 136
gribcodes_index_select(iiduvidx, 'levelshortName',lev param) #136
gid grib_index_select(iiduv,'shortName',param= codes_new_from_index(idx)
values[param] = codes_get_values(gid)
gidp=grib_new_from_index(iiduv codes_release(gid)
if lev values_lev[param] = grib_get_values(gidp)< values['nlevels']:
grib_release(gidp)
codes_index_select(idx, 'level', lev + 1) # 137
gid if (lev < levelSize):= codes_new_from_index(idx)
values['prev' + param] grib_index_select(iiduv,'level',lev+1) #137= codes_get_values(gid)
gidp=grib_new_from_index(iiduvcodes_release(gid)
else:
values_plev[['prev' + param] = grib_get_values(gidpnp.zeros(values['sp'].size)
# search if the provided wind height converted to
# geopotential (my_z) grib_release(gidpis between the current level (z_f)
# and the previous one (z_f_prev)
for i elsein range(z_f_prev.size):
if found[i]:
values_plev[param] = zeros(sp.size) continue
for i in arange(z_f_prev.size)if values['gh'][i] >= z_f_prev[i] and values['gh'][i] < z_f[i]:
found[i] = True
if found[i]: continue
# when found, interpolate vertically to get the
if my_z[i] >= z_f_prev[i] and my_z[i] < z_f[i]:
# value and store it in out[param] to be written
# at the end
found[i]=True
for param in params:
#print "wind %d for point %d found betweenres ml %d(%lf) and %d(%lf)\n" %(my_z[i],i,lev+1,z_f_prev[i],lev,z_f[i])
= (((float(values[param][i]) *
(values['gh'][i] for param in ["u","v"]:- z_f_prev[i])) +
(float(values['prev' + out[paramparam][i] = ((values_lev[param][i] * ( my_z[i]-z_f_prev[i])) + (values_plev[param][i] *) *
(z_f[i] - my_zvalues['gh'][i]) )) /
(z_f[i] - z_f_prev[i]))
z_f_prev=z_f out[param][i] = res
# update z_f_prev
for i in arange(sp.size): z_f_prev = z_f
# simple error check
for i in range(values['sp'].size):
if not found[i]:
print "point ",i,"print('point ', i, 'not found..."')
# write the values in the fout file
for param in ["u","v"]params:
codes_set(values['sample'], 'shortName', param)
gribcodes_set_(values(gid_out[param],out[param]['sample'], 'typeOfLevel', 'heightAboveGround')
codes_set(values['sample'], 'level', values['height'])
grib_write(gid_codes_set_values(values['sample'], out[param],fout)
grib_release(gid_out[param])
grib_index_release(iidtq)
grib_index_release(iiduv)
grib_index_release(zlnsp)
fout.close()
print("Done")
if __name__ == "__main__":
request_date=0
request_time=0
wind = 100
parser = argparse.ArgumentParser(
description='Python tool to calculate the Z of the model levels')
parser.add_argument("-w","--wind", help="height to calculate the wind components",required=True)
parser.add_argument("-o","--output", help="name of the output file")
parser.add_argument('t_q', metavar='tq.grib', type=str,
help='grib file with temperature(t) and humidity(q) for the model levels')
parser.add_argument('z_lnsp', metavar='zlnsp.grib', type=str,
help='grib file with geopotential(z) and Logarithm of surface pressure(lnsp) for the ml=1')
parser.add_argument('u_v', metavar='uv.grib', type=str,
help='grib file with u and v component of wind for the model levels')
args = parser.parse_args()
for fname in (args.t_q,args.z_lnsp,args.u_v):
if not os.path.isfile(fname):
print "[ERROR] file %s does not exist" %(fname)
sys.exit(1)
if args.wind:
wind = args.wind
out_name = 'uv_out_'+wind+'m.grib'
if args.output:
out_name=args.output
#calling main function
main(args.t_q,args.z_lnsp,args.u_v,out_name,wind)
codes_write(values['sample'], fout)
class WrongStepError(Exception):
''' Exception capturing wrong step'''
pass
if __name__ == '__main__':
main()
|