Description

This example shows: how to set pv values in a GRIB message

Source code

grib_set_pv.f90
! (C) Copyright 2005- 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.
!
!  Description: how to set pv values in a GRIB message
!
!
program grib_set_pv
  use eccodes
  implicit none
  integer                         :: numberOfLevels
  integer                         :: numberOfCoefficients
  integer                         :: outfile, igrib
  integer                         :: i, ios
  real, dimension(:),allocatable  :: pv

  numberOfLevels=60
  numberOfCoefficients=2*(numberOfLevels+1)

  allocate(pv(numberOfCoefficients))

  ! read the model level coefficients from file
  open( unit=1, file="../../data/60_model_levels", &
                form="formatted",action="read")

  do i=1,numberOfCoefficients,2
    read(unit=1,fmt=*, iostat=ios) pv(i), pv(i+1)
    if (ios /= 0) then
      print *, "I/O error: ",ios
      exit
    end if
  end do

  ! print coefficients
  !do i=1,numberOfCoefficients,2
  !  print *,"  a=",pv(i)," b=",pv(i+1)
  !end do

  close(unit=1)

  call codes_open_file(outfile, 'out.pv.grib1','w')

  ! A new grib message is loaded from file
  ! igrib is the grib id to be used in subsequent calls
  call codes_grib_new_from_samples(igrib, "reduced_gg_sfc_grib1")

  ! set levtype to ml (model level)
  call codes_set(igrib,'typeOfLevel','hybrid')

  ! set level
  call codes_set(igrib,'level',2)

  ! set PVPresent as an integer
  call codes_set(igrib,'PVPresent',1)

  call codes_set(igrib,'pv',pv)

  ! write modified message to a file
  call codes_write(igrib,outfile)

  ! free memory
  call codes_release(igrib)
  deallocate(pv)

  call codes_close_file(outfile)

end program grib_set_pv
grib_set_pv.py
#
# (C) Copyright 2005- 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.
#
# Description: how to set pv values in a GRIB message
#

import traceback
import sys

from eccodes import *

VERBOSE = 1  # verbose error reporting


def example():
    """
    Encoding of the pv coefficients.
    """
    # read the coefficients from file
    pv = []
    for line in open('../../data/60_model_levels'):
        pv.extend([float(x) for x in line.strip().split('\t')])

    numberOfLevels = 60
    numberOfCoefficients = 2 * (numberOfLevels + 1)
    assert (len(pv) == numberOfCoefficients)

    fout = open('grib_set_pv.py.temp.grib', 'wb')
    gid = codes_grib_new_from_samples('reduced_gg_sfc_grib1')

    codes_set(gid, 'typeOfLevel', 'hybrid')
    codes_set(gid, 'level', 2)
    codes_set(gid, 'PVPresent', 1)
    codes_set_array(gid, 'pv', pv)

    codes_write(gid, fout)

    codes_release(gid)
    fout.close()


def main():
    try:
        example()
    except CodesInternalError as err:
        if VERBOSE:
            traceback.print_exc(file=sys.stderr)
        else:
            sys.stderr.write(err.msg + '\n')

        return 1


if __name__ == '__main__':
    sys.exit(main())
grib_set_pv.c
/*
 * (C) Copyright 2005- 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.
 */

/*
 * C Implementation: grib_set_pv
 *
 * Description: how to set pv (vertical coordinate parameters)
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "eccodes.h"

static void usage(const char* prog) {
    fprintf(stderr, "usage: %s in out\n",prog);
    exit(1);
}

int main(int argc, char** argv)
{
    int err = 0;
    long NV = 0;
    size_t size=0;
    double pv[4]={1,2,3,4};
    size_t pvsize=4;

    FILE* in = NULL;
    char* infile = NULL;
    FILE* out = NULL;
    char* outfile = NULL;
    codes_handle *h = NULL;
    const void* buffer = NULL;

    if (argc != 3) usage(argv[0]);
    infile = argv[1];
    outfile= argv[2];

    in = fopen(infile, "r");
    if(!in) {
        fprintf(stderr, "ERROR: unable to open input file %s\n",infile);
        return 1;
    }

    out = fopen(outfile, "w");
    if(!out) {
        fprintf(stderr, "ERROR: unable to open output file %s\n",outfile);
        fclose(in);
        return 1;
    }

    /* create a new handle from a message in a file */
    h = codes_handle_new_from_file(0, in, PRODUCT_GRIB, &err);
    if (h == NULL) {
        fprintf(stderr, "Error: unable to create handle from file %s\n",infile);
    }

    CODES_CHECK(codes_set_long(h,"PVPresent", 1),0);

    CODES_CHECK(codes_set_double_array(h, "pv", pv, pvsize),0);

    /* Once we set the pv array, the NV key should be also set */
    CODES_CHECK(codes_get_long(h, "NV", &NV),0);
    assert( NV == pvsize );

    /* get the coded message in a buffer */
    CODES_CHECK(codes_get_message(h, &buffer, &size),0);

    /* write the buffer in a file*/
    if(fwrite(buffer, 1, size, out) != size)
    {
        perror(outfile);
        exit(1);
    }

    codes_handle_delete(h);
    fclose(in);
    fclose(out);

    return 0;
}