bufr_read_temp

! (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.
!
!
! FORTRAN 90 Implementation: bufr_read_temp
!
! Description: how to read levels from TEMP BUFR messages.
!
! Please note that TEMP reports can be encoded in various ways in BUFR. Therefore the code
! below might not work directly for other types of TEMP messages than the one used in the
! example. It is advised to use bufr_dump first to understand the structure of these messages.
!
program bufr_read_temp
  use eccodes
  implicit none
  integer            :: ifile
  integer            :: iret
  integer            :: ibufr
  integer            :: i, count=0
  integer(kind=4),dimension(:), allocatable  :: timePeriod,extendedVerticalSoundingSignificance
  integer(kind=4)  :: blockNumber,stationNumber
  real(kind=8),dimension(:), allocatable :: pressure,airTemperature,dewpointTemperature
  real(kind=8),dimension(:), allocatable :: geopotentialHeight,latitudeDisplacement,longitudeDisplacement
  real(kind=8),dimension(:), allocatable :: windDirection,windSpeed

  call codes_open_file(ifile,'../../data/bufr/PraticaTemp.bufr','r')

  ! The first bufr message is loaded from file,
  ! ibufr is the bufr id to be used in subsequent calls
  call codes_bufr_new_from_file(ifile,ibufr,iret)
  do while (iret/=CODES_END_OF_FILE)
    write(*,*) 'message: ',count
    call codes_set(ibufr,'unpack',1)
    call codes_get(ibufr,'timePeriod',timePeriod)
    call codes_get(ibufr,'pressure',pressure)
    call codes_get(ibufr,'extendedVerticalSoundingSignificance',extendedVerticalSoundingSignificance)
    call codes_get(ibufr,'nonCoordinateGeopotentialHeight',geopotentialHeight)
    call codes_get(ibufr,'latitudeDisplacement',latitudeDisplacement)
    call codes_get(ibufr,'longitudeDisplacement',longitudeDisplacement)
    call codes_get(ibufr,'airTemperature',airTemperature)
    call codes_get(ibufr,'dewpointTemperature',dewpointTemperature)
    call codes_get(ibufr,'windDirection',windDirection)
    call codes_get(ibufr,'windSpeed',windSpeed)
    call codes_get(ibufr,'blockNumber',blockNumber)
    call codes_get(ibufr,'stationNumber',stationNumber)
    print *,'station',blockNumber,stationNumber
    print *,'timePeriod pressure geopotentialHeight latitudeDisplacement &
          &longitudeDisplacement airTemperature windDirection windSpeed significance'
    do i=1,size(windSpeed)
      write(*,'(I5,6X,F9.1,2X,F9.2,10X,F8.2,14X,F8.2,16X,F8.2,6X,F8.2,4X,F8.2,4X,I0)') timePeriod(i),pressure(i),&
          &geopotentialHeight(i),latitudeDisplacement(i),&
          &longitudeDisplacement(i),airTemperature(i),windDirection(i),windSpeed(i),extendedVerticalSoundingSignificance(i)
    enddo
    ! Free arrays
    deallocate(timePeriod)
    deallocate(pressure)
    deallocate(geopotentialHeight)
    deallocate(latitudeDisplacement)
    deallocate(longitudeDisplacement)
    deallocate(airTemperature)
    deallocate(dewpointTemperature)
    deallocate(windDirection)
    deallocate(windSpeed)
    deallocate(extendedVerticalSoundingSignificance)
    ! Release the bufr message
    call codes_release(ibufr)
    ! Load the next bufr message
    call codes_bufr_new_from_file(ifile,ibufr,iret)
    count=count+1
  end do
  ! Close file
  call codes_close_file(ifile)

end program bufr_read_temp

#
# (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.
#
# Python implementation: bufr_read_temp
#
# Description: how to read levels from TEMP BUFR messages.
#
#
#
# Please note that TEMP reports can be encoded in various ways in BUFR.
# Therefore the code below might not work directly for other types of TEMP
# messages than the one used in the example. It is advised to use bufr_dump to
# understand the structure of the messages.
#
import sys
import traceback

from eccodes import *

INPUT = '../../data/bufr/PraticaTemp.bufr'
VERBOSE = 1  # verbose error reporting


def example():
    # open BUFR file
    f = open(INPUT, 'rb')
    cnt = 0
    # loop over the messages in the file
    while 1:
        # get handle for message
        bufr = codes_bufr_new_from_file(f)
        if bufr is None:
            break
        print("message: %s" % cnt)
        # we need to instruct ecCodes to expand all the descriptors
        # i.e. unpack the data section
        codes_set(bufr, 'unpack', 1)
        # get all the timePeriods
        timePeriod = codes_get_array(bufr, "timePeriod")
        pressure = codes_get_array(bufr, "pressure")
        extendedVerticalSoundingSignificance = codes_get_array(bufr, "extendedVerticalSoundingSignificance")
        geopotentialHeight = codes_get_array(bufr, "nonCoordinateGeopotentialHeight")
        latitudeDisplacement = codes_get_array(bufr, "latitudeDisplacement")
        longitudeDisplacement = codes_get_array(bufr, "longitudeDisplacement")
        airTemperature = codes_get_array(bufr, "airTemperature")
        dewpointTemperature = codes_get_array(bufr, "dewpointTemperature")
        windDirection = codes_get_array(bufr, "windDirection")
        windSpeed = codes_get_array(bufr, "windSpeed")
        blockNumber = codes_get(bufr, "blockNumber")
        stationNumber = codes_get(bufr, "stationNumber")
        print('station %d%d' % (blockNumber, stationNumber))
        print(
            'timePeriod pressure geopotentialHeight latitudeDisplacement longitudeDisplacement airTemperature windDirection windSpeed significance')
        for i in range(0, len(windSpeed) - 1):
            print(timePeriod[i], pressure[i], geopotentialHeight[i], latitudeDisplacement[i], longitudeDisplacement[i],
                  airTemperature[i], windDirection[i], windSpeed[i], extendedVerticalSoundingSignificance[i])
        cnt += 1
        # delete handle
        codes_release(bufr)
    # close the file
    f.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())

/*
 * (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: bufr_read_temp
 *
 * Description: how to read temperature significant levels from TEMP BUFR messages.
 *
 */

/*
 * Please note that TEMP reports can be encoded in various ways in BUFR. Therefore the code
 * below might not work directly for other types of TEMP messages than the one used in the
 * example. It is advised to use bufr_dump to understand the structure of the messages.
 */

#include "eccodes.h"

int main(int argc,char* argv[])
{
    FILE* in = NULL;

    /* Message handle. Required in all the eccodes calls acting on a message.*/
    codes_handle* h=NULL;

    double *sigt_pres=NULL, *sigt_geo=NULL, *sigt_t=NULL;
    double *sigt_td=NULL;
    size_t sigt_len=0, len=0, i=0;
    int err=0;
    int cnt=0;
    const char* infile = "../../data/bufr/temp_101.bufr";
    char key_name[128];

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

    /* Loop over the messages in the bufr file */
    while ((h = codes_handle_new_from_file(NULL,in,PRODUCT_BUFR,&err)) != NULL || err != CODES_SUCCESS)
    {
        if (h == NULL) {
            printf("Error: unable to create handle for message %d\n",cnt);
            cnt++;
            continue;
        }

        printf("message: %d\n",cnt);

        /* We need to instruct ecCodes to expand the descriptors
           i.e. unpack the data values */
        CODES_CHECK(codes_set_long(h,"unpack",1),0);

        /* In what follows we rely on the fact that for
         * temperature significant levels the value of key
         * verticalSoundingSignificance is 4 (see flag table 8001 for details).
         *
         * In our BUFR message verticalSoundingSignificance is always followed by
         *    geopotential, airTemperature, dewpointTemperature,
         *    windDirection, windSpeed and pressure.
         * So in order to access any of these keys we need to use the
         * condition: verticalSoundingSignificance=4.
         */

        /* Get the number of the temperature significant levels.*/

        /* We find out the number of temperature significant levels by
         * counting how many pressure values we have on these levels.*/

        sprintf(key_name,"/verticalSoundingSignificance=4/pressure");
        CODES_CHECK(codes_get_size(h,key_name,&sigt_len),0);

        printf("Number of T significant levels: %lu\n",sigt_len);

        /* Allocate memory for the values to be read. Each
         * parameter must have the same number of values. */
        sigt_pres = (double*)malloc(sigt_len*sizeof(double));
        sigt_geo = (double*)malloc(sigt_len*sizeof(double));
        sigt_t = (double*)malloc(sigt_len*sizeof(double));
        sigt_td = (double*)malloc(sigt_len*sizeof(double));

        /* Get pressure */
        sprintf(key_name,"/verticalSoundingSignificance=4/pressure");
        len=sigt_len;
        CODES_CHECK(codes_get_double_array(h,key_name,sigt_pres,&len),0);

        /* Get geopotential */
        sprintf(key_name,"/verticalSoundingSignificance=4/nonCoordinateGeopotential");

        /* Check the size */
        CODES_CHECK(codes_get_size(h,key_name,&len),0);
        if(len != sigt_len)
        {
            printf("inconsistent number of geopotential values found!\n");
            return 1;
        }

        /* Get the values */
        CODES_CHECK(codes_get_double_array(h,key_name,sigt_geo,&len),0);

        /* Get temperature */
        if(len != sigt_len) /* Check the size */
        {
            printf("inconsistent number of temperature values found!\n");
            return 1;
        }

        /* Get the values */
        sprintf(key_name,"/verticalSoundingSignificance=4/airTemperature");
        CODES_CHECK(codes_get_double_array(h,key_name,sigt_t,&len),0);

        /* Get dew point */
        if(len != sigt_len) /* Check the size */
        {
            printf("inconsistent number of dewpoint temperature values found!\n");
            return 1;
        }

        /* Get the values */
        sprintf(key_name,"/verticalSoundingSignificance=4/dewpointTemperature");
        CODES_CHECK(codes_get_double_array(h,key_name,sigt_td,&len),0);

        /* Print the values */
        printf("lev  pres    geo    t    td\n");
        printf("-------------------------------\n");

        for(i=0; i < sigt_len; i++)
        {
            printf("%3lu %6.0f %6.0f %.1f %.1f\n",
                    i+1,sigt_pres[i],sigt_geo[i],sigt_t[i],sigt_td[i]);
        }

        /* Delete handle */
        codes_handle_delete(h);

        /* Release memory */
        free(sigt_pres);
        free(sigt_geo);
        free(sigt_t);
        free(sigt_td);

        cnt++;
    }

    fclose(in);
    return 0;
}