Versions Compared

Old Version 4

changes.mady.by.user Shahram Najm

Saved on

compared with

New Version 5

changes.mady.by.user Xavier Abellan

Saved on

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.

...

Code Block
languagepy
titlecompute_wind_speed_height.py
#!/usr/bin/env python
#'''
# Copyright 20152018 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_wind_speed_height
#
# Author (date) : Cristian Simarro (20/11/2015)
#
# modified:       Xavi Abellan     (03/12/2018) - compatibilty with Python 3

Category      : COMPUTATION
#
# OneLineDesc   : Computes the u/v components of the wind at certain height
#
# Description   : Computes computes the u/v components of the wind at certain height.
#                height. First, it calculates the geopotential of each model
                level. Once the requested height 
#is found between two model
             is found between two model levels, the program will vertically interpolate the u/v
                component values.
#                 Based on code from Nils Wedi, the IFS documentation:
#https://www.ecmwf.int/en/forecasts/documentation-and-support/changes-ecmwf-model/ifs-documentation
                part https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation
#                 part III. III. Dynamics and numerical procedures
#                 optimised implementation by Dominique Lucas.
#                 ported to Python by Cristian Simarro
#
# Parameters    : -w wind                - height in meters you want to know u/v components
#                 tq.grib                - grib file with all the levelist of t and q
#u/v components
                 uvtq.grib                - grib file with all the levelistslevelist of u/v
#
                    zlnsp.grib             - grib file with levelist 1 for params zt and lnspq
#                uv.grib -o output   (optional) - name of the output file (default='uv_out_<wind>.grib')
#
# Return Value  : output (default='uv_out_<wind>.grib')
#      - grib file with all the levelists of
           A fieldset the u/v components values for at the specified height
#
# Dependencies  : None
#
# Example Usage : 
#                                  u/v
             python compute_wind_speed_height.py tq_ml_20151116_0.grib zlnsp_ml_20151116_0.grib uv_ml_20151116_0.grib -w 100

from numpy import * 
import sys,math,os
import argparse 
from gribapi import *

PARAMS = ["u","v"]

def main(t_q,z_lnsp,u_v,out_name,h):

       - grib file with levelist 1 for
      #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):
     params z and lnsp
    os.remove(out_name)
    fout = open(out_name,'w')
    ftmp = open(t_q)
  -o output   ftmp.close(optional)
 - name of Rdthe =output 287.06file
     RG = 9.80665
    index_keys = ["date","time","shortName","level","step"]

    values= {}
    pv = {}
    out = {}
    gid_out = {}
    values_plev = {}
    values_lev = {}
(default='uv_out_<wind>.grib')

Return Value  : output (default='uv_out_<wind>.grib')
       
    h=int(h)
    zlnsp = grib_index_new_from_file(z_lnsp,index_keys)
    iidtq = grib_index_new_from_file(t_q,index_keys)
A fieldset the u/v components values for at the specified
      iiduv = grib_index_new_from_file(u_v,index_keys)
    
    counter=0height

Dependencies  : None

Example #Usage iterate:
 date
    for date in grib_index_get(zlnsp,'date'):
        gribcompute_wind_index_select(zlnsp,'date',date)
        grib_index_select(iidtq,'date',date)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,
        grib             codes_index_select(iiduv,'date',date), codes_new_from_index, codes_set,
        # iterate step
        for time in gribcodes_index_get(zlnsp,'time'):add_file, codes_get_array, codes_get_values,
            grib_index_select(zlnsp,'time',time)
            grib_index_select(iidtq,'time',time)codes_index_release, codes_release, codes_set_values,
            grib_index_select(iiduv,'time',time)
         codes_write)

R_D =  grib_index_select(zlnsp,'level',1)287.06
R_G = 9.80665


def parse_args():
    ''' Parse program arguments using ArgumentParser'''
   grib_index_select(zlnsp,'step',0) parser = argparse.ArgumentParser(
        description='Python tool to calculate the  grib_index_select(zlnsp,'shortName','z')wind at certain height')
    parser.add_argument('-w', '--height', required=True, type=int,
            gid = grib_new_from_index(zlnsp)

          help='height to #gridTypecalculate mustthe be gridded, not spectralwind components')
            if grib_get(gid,"gridType",str) == "sh":
         parser.add_argument('-o', '--output', help='name of the output file')
    parser.add_argument('t_q', metavar='tq.grib', type=str,
       print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral')
          help=('grib file with temperature(t) and  sys.exit(1)

humidity(q)',
              # surface geopotential
            values["z"] = grib_get_values(gid)
     'for the model levels'))
    parser.add_argument('z_lnsp', metavar='zlnsp.grib', type=str,
        z_h  = values["z"]
            pv help= grib_get_array(gid,'pv')('grib file with geopotential(z) and Logarithm',
            levelSizeNV = grib_get(gid,'NV',int)/2 -1
            
   'of surface pressure(lnsp) for the ml=1'))
    gribparser.add_release(gid)argument('u_v', metavar='uv.grib', type=str,
           
            # iterate step all but geopotential z which is always step 0 (an)help=('grib file with u and v component of wind for',
            for step in grib_index_get(iidtq,'step'):
               'the # we need to get z and lnsp from the first level to do the calculationsmodel levels'))
    args = parser.parse_args()
    if not args.output:
        args.output        z_h = values["z"]= 'uv_out_%sm.grib' % args.height
    return args


def main():
    '''Main function'''
    args # heigh in geopotential
      = parse_args()

    print('Arguments: %s' % ", ".join(
        ['%s: %s' my_z = h*RG + z_h
   % (k, v) for k, v in vars(args).items()]))

    fout = open(args.output, 'wb')
      z_f_previndex_keys = z_h
    ['date', 'time', 'shortName', 'level', 'step']

    idx        = codes_index_new_from_file(args.z_lnsp, index_keys)
    codes_index_add_file(idx, args.t_q)
    if 'u_v' in args:
         gribcodes_index_selectadd_file(iidtq,'step',stepidx, args.u_v)

    # iterate date
    for date in    gribcodes_index_selectget(iiduvidx, 'stepdate',step):

        codes_index_select(idx, 'date', date)
       for shortName# in ["lnsp"]:iterate step
        for time          gribin codes_index_selectget(zlnspidx, 'shortNametime',shortName):
                   gribcodes_index_select(zlnspidx, 'steptime',step time)
            values       gid = grib_new_from_index(zlnsp)
= get_initial_values(idx, keep_sample=True)
            if 'height' in args:
        if grib_get(gid,"gridType",str) == "sh":
     values['height'] = args.height
                print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral')values['gh'] = args.height * R_G + values['z']
            if 'levelist' in args:
        sys.exit(1)
        values['levelist'] = args.levelist
           values[shortName] = grib_get_values(gid)
          # iterate step all but geopotential z which is always step 0 (an)
          pv = grib_get_array(gid,'pv') for step in codes_index_get(idx, 'step'):
                   levelSizeNV = grib_get(gid,'NV',int)/2 -1codes_index_select(idx, 'step', step)
                #   grib_release(gid)surface pressure
   
             values['sp'] =  # surface pressureget_surface_pressure(idx)
                sp = exp(values["lnsp"]production_step(idx, values, fout)

        try:
        # get the coefficients for computing the pressures codes_release(values['sample'])
        except KeyError:
       # how many levels are we computing? pass

                gribcodes_index_select(iidtq,'shortName',"t"release(idx)
                levelSize=max(grib_index_get(iidtq,"level",int))fout.close()


def get_initial_values(idx, keep_sample=False):
    '''Get the values of surface z, pv and number of   if levelSize != levelSizeNV:levels '''
    codes_index_select(idx, 'level', 1)
        codes_index_select(idx, 'step', 0)
    codes_index_select(idx,  print(sys.argv[0]+' [WARN] total levels should be: '+str(levelSizeNV)+' but it is '+str(levelSize))
                A = pv[0:levelSize+1]
                B = pv[levelSize+1:]
                Ph_levplusone = A[levelSize] + (B[levelSize]*sp)

    'shortName', 'z')
    gid = codes_new_from_index(idx)

    values = {}
    # surface geopotential
    values['z'] = codes_get_values(gid)
    values['pv'] = codes_get_array(gid, 'pv')
    values['nlevels'] = codes_get(gid, 'NV', int) // 2 - 1
    check_max_level(idx, values)
    if keep_sample:
        values['sample'] = gid
    else:
        codes_release(gid)
    # We want to integrate up into the atmosphere, starting at the groundreturn values


def check_max_level(idx, values):
    '''Make sure we have all the levels required'''
    # how many levels are we computing?
    max_level  # so we start at the lowest level (highest number) and keep= max(codes_index_get(idx, 'level', int))
    if max_level != values['nlevels']:
        print('%s [WARN] total levels should be: %d but #it accumulatingis the height as we go.%d' %
              (sys.argv[0],  # See the IFS documentation:values['nlevels'], max_level),
                # https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentationfile=sys.stderr)
             values['nlevels'] = max_level


def get_surface_pressure(idx):
   # part III
    '''Get the surface pressure for date-time-step'''
    codes_index_select(idx, 'level', 1)
    codes_index_select(idx, 'shortName', 'lnsp')
 # For speed andgid file I/O, we perform the computations with numpy vectors instead= codes_new_from_index(idx)
    if codes_get(gid, 'gridType', str) == 'sh':
        print('%s [ERROR] fields must be gridded, not spectral' #% of fieldsets.sys.argv[0],
                file=sys.stderr)
        sys.exit(1)
    # surface pressure
  #initialize values forsfc_p the output= np.exp(codes_get_values(gid))
    codes_release(gid)
    return sfc_p


def       for param in PARAMSget_ph_levs(values, level):
    '''Return the presure at a given level and the next'''
    a_coef   grib_index_select(iiduv,'level',1)= values['pv'][0:values['nlevels'] + 1]
    b_coef = values['pv'][values['nlevels'] + 1:]
    ph_lev = a_coef[level - 1] + (b_coef[level - 1]  grib_index_select(iiduv,'shortName',param* values['sp'])
    ph_levplusone = a_coef[level] + (b_coef[level] * values['sp'])
    return ph_lev, ph_levplusone


def compute_z_level(idx, lev,  gid = grib_new_from_index(iiduv)values, z_h):
    '''Compute z at half & full level for the given level, based on t/q/sp'''
   gid_out[param]=grib_clone(gid)
             # select the levelist and retrieve the vaules of t and q
    # t_level: values  grib_release(gid)for t
    # q_level: values for q
    codes_index_select(idx, 'level', lev)
    codes_index_select(idx,  out[param]=zeros(sp.size'shortName', 't')
    gid = codes_new_from_index(idx)
          found = [False for i in range(sp.size)]t_level = codes_get_values(gid)
    codes_release(gid)
    codes_index_select(idx, 'shortName', 'q')
    gid = codes_new_from_index(idx)
    forq_level lev in list(reversed(range(1,levelSize+1))):= codes_get_values(gid)
    codes_release(gid)

    # compute moist temperature
    t_level = t_level * (1. #+ select0.609133 the levelist and retrieve the vaules of t and q* q_level)

    # compute the pressures (on half-levels)
    ph_lev, ph_levplusone = get_ph_levs(values, lev)

    if lev == 1:
     #   t_level: values for tdlog_p = np.log(ph_levplusone / 0.1)
        alpha = np.log(2)
    else:
       # q_level: values for qdlog_p = np.log(ph_levplusone / ph_lev)
        alpha = 1. - ((ph_lev / (ph_levplusone - ph_lev)) *   grib_index_select(iidtq,'level',levdlog_p)

    t_level = t_level * R_D

    # z_f is the geopotential of this  grib_index_select(iidtq,'shortName',"t")full level
    # integrate from previous (lower) half-level z_h to the
    # full level
    gidz_f = grib_new_from_index(iidtq)z_h + (t_level * alpha)

    # z_h is the geopotential of 'half-levels'
    # integrate z_h to next half t_level = grib_get_values(gid)
    z_h = z_h + (t_level * dlog_p)

    return z_h, z_f


def production_step(idx,   grib_release(gid)values, fout):
    '''Produce u/v interpolated from ML for a given  height'''
    # We  grib_index_select(iidtq,'shortName',"q")
        want to integrate up into the atmosphere, starting at the
    # ground so we start at the lowest gidlevel = grib_new_from_index(iidtq)(highest number) and
    # keep accumulating the height as we go.
    # See the IFS documentation, q_level = grib_get_values(gid)
part III
    # For speed and file I/O, we perform the computations with
    # numpy vectors instead  grib_release(gid)of fieldsets.

    out = {}
    params =     ('u', 'v')
    # computewe moistneed temperature
to get z and lnsp from the first level to do the
    # calculations
    tz_levelh = t_level * (1.+0.609133*q_level)

values['z']
    # height in geopotential
    z_f_prev = z_h
    # initialize values #for compute the pressures (on half-levels)output
    for param in params:
             Ph_lev = A[lev-1] + (B[lev-1] * sp)

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'] + if lev == 11)))):
        z_h, z_f = compute_z_level(idx, lev,            dlogP = log(Ph_levplusone/0.1values, z_h)
        # retrieve u/v params for the current level
        for alphaparam = log(2)in params:
            codes_index_select(idx, 'level', lev)  #    else:136
                codes_index_select(idx, 'shortName', param)
            dlogPgid = log(Ph_levplusone/Ph_levcodes_new_from_index(idx)
            values[param] = codes_get_values(gid)
          dP  codes_release(gid)
  = Ph_levplusone-Ph_lev
         if lev < values['nlevels']:
            alpha = 1. - ((Ph_lev/dP)*dlogP)

   codes_index_select(idx, 'level', lev + 1)  # 137
                gid TRd = t_level*Rd

codes_new_from_index(idx)
                values['prev' + param]  # z_f is the geopotential of this full level
= codes_get_values(gid)
                codes_release(gid)
           # integrateelse:
 from previous (lower) half-level z_h to the full level
       values['prev' + param] = np.zeros(values['sp'].size)
        # z_fsearch =if z_hthe + (TRd*alpha)

   provided wind height converted to
        # geopotential (my_z) is between the current level  # (z_hf)
 is the geopotential of 'half-levels'
   # and the previous one (z_f_prev)
        for i   # integrate z_h to next half levelin range(z_f_prev.size):
            if found[i]:
       z_h=z_h+(TRd*dlogP)

         continue
           Ph_levplusone = Ph_lev
      if values['gh'][i] >= z_f_prev[i] and values['gh'][i] < z_f[i]:
              # retrieve u/v params for the current level found[i] = True
                # when found, interpolate forvertically paramto inget PARAMS:the
                # value and store it in out[param]  grib_index_select(iiduv,'level',lev) #136
to be written
                #  at the end
         grib_index_select(iiduv,'shortName',param)
         for param in params:
                    res  gidp=grib_new_from_index(iiduv)= (((float(values[param][i]) *
                        values_lev[param] = grib_get_values(gidp)     (values['gh'][i] - z_f_prev[i])) +
                        grib_release(gidp)
                        if (lev < levelSize):
  (float(values['prev' + param][i]) *
                          grib_index_select(iiduv,'level',lev+1) #137
                            gidp=grib_new_from_index(iiduv)
                            values_plev[param] = grib_get_values(gidp)
                            grib_release(gidp)
                        else:
                            values_plev[param] = zeros(sp.size)
                    # search if the provided wind height converted to geopotential (my_z) is between
                    # the current level (z_f) and the previous one (z_f_prev)
                    for i in arange(z_f_prev.size):
                        if found[i]: continue
                        if my_z[i] >= z_f_prev[i] and my_z[i] < z_f[i]:
                            found[i]=True
                            # when found, interpolate vertically to get the value
                            # store it in out[param] to be written at the end
                            for param in PARAMS:
                                out[param][i] = ( (values_lev[param][i] * ( my_z[i]-z_f_prev[i])) + \
                                                 (values_plev[param][i] * (z_f[i] - my_z[i]) ) ) /  \
                                                 (z_f[i] - z_f_prev[i])
                    # update z_f_prev
                    z_f_prev=z_f

                # simple error check
                for i in arange(sp.size):
                    if not found[i]:
                        print "point ",i,"not found..."

                # write the values in the fout file
                for param in PARAMS:
                    grib_set_values(gid_out[param],out[param])
                    grib_write(gid_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, 
            (z_f[i] - values['gh'][i]))) /
              help='grib file with temperature(t) and humidity(q) for the model levels')
    parser.add_argument('z_lnsp', metavar='zlnsp.grib', type=str, f[i] - z_f_prev[i]))
                   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, 
      out[param][i] = res
        # update z_f_prev
        z_f_prev = z_f

    help='grib file with u and v component of wind for the model levels')# simple error check
    for i in range(values['sp'].size):
    args = parser.parse_args()
  if  for fname in (args.t_q,args.z_lnsp,args.u_v)not found[i]:
        if not os.path.isfile(fname):    print('point ', i, 'not found...')

    # write the values in the fout file
   print "[ERROR] file %s does not exist" %(fname)
  for param in params:
        codes_set(values['sample'], 'shortName', param)
        codes_set(values['sample'],   sys.exit(1'typeOfLevel', 'heightAboveGround')
    if args.wind:
   codes_set(values['sample'], 'level', values['height'])
   wind = args.wind
    out_name = 'uv_out_'+wind+'m.grib'
    if args.output:codes_set_values(values['sample'], out[param])
        out_name=args.output

    #calling main functioncodes_write(values['sample'], fout)


if __name__ == '__main__':
    main(args.t_q,args.z_lnsp,args.u_v,out_name,wind))