Page History

!Copyright 2005-2018 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_tropical_cyclone
!
! Description: how to read data for a tropical cyclone BUFR message.
!
! Please note that tropical cyclone tracks can be encoded in various ways in BUFR.
! Therefore the code below might not work directly for other types of messages
! than the one used in the example. It is advised to use bufr_dump to
! understand the structure of the messages.
!
program bufr_read_tropical_cyclone
  use eccodes
  implicit none
  integer            :: ifile
  integer            :: iret
  integer            :: ibufr,skipMember
  integer            :: significance
  integer            :: year,month,day,hour,minute
  integer            :: i,j,k,ierr,count=1
  integer            :: rankPosition,rankSignificance,rankPressure,rankWind
  integer            :: rankPeriod,numberOfPeriods

  real(kind=8) :: latitudeCentre,longitudeCentre
  real(kind=8), dimension(:), allocatable :: latitudeMaxWind0,longitudeMaxWind0,windMaxWind0
  real(kind=8), dimension(:), allocatable :: latitudeAnalysis,longitudeAnalysis,pressureAnalysis
  real(kind=8), dimension(:,:), allocatable :: latitude,longitude,pressure
  real(kind=8), dimension(:,:), allocatable :: latitudeWind,longitudeWind,wind
  integer(kind=4), dimension(:), allocatable :: memberNumber,period
  real(kind=8), dimension(:), allocatable :: values
  integer(kind=4), dimension(:), allocatable :: ivalues

  character(len=8)   :: rankSignificanceStr,rankPositionStr,rankPressureStr,rankWindStr
  character(len=8)   :: stormIdentifier,rankPeriodStr

  call codes_open_file(ifile,'../../data/bufr/tropical_cyclone.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(*,'(A,I3,A)') '**************** MESSAGE: ',count,'  *****************'

    ! We need to instruct ecCodes to unpack the data values
    call codes_set(ibufr,"unpack",1);

    call codes_get(ibufr,'year',year);
    call codes_get(ibufr,'month',month);
    call codes_get(ibufr,'day',day);
    call codes_get(ibufr,'hour',hour);
    call codes_get(ibufr,'minute',minute);
    write(*,'(A,I0,A,I0,A,I0,A,I0,A,I0,A,I0)')'Date and time: ',day,'.',month,'.',year,'  ',hour,':',minute

    call codes_get(ibufr,'stormIdentifier',stormIdentifier)
    write(*,'(A,A)')'Storm identifier: ',stormIdentifier

    ! How many different timePeriod in the data structure?
    rankPeriod=0
    ierr=0
    do while(ierr==0)
      rankPeriod=rankPeriod+1
      write (rankPeriodStr,'(I0)')rankPeriod
      call codes_get(ibufr,'#'//trim(rankPeriodStr)//'#timePeriod',period,ierr)
      if(allocated(period)) deallocate(period)
    enddo
    ! The numberOfPeriods includes the analysis (period=0)
    numberOfPeriods=rankPeriod

    call codes_get(ibufr,'ensembleMemberNumber',memberNumber)

    allocate(latitude(size(memberNumber),numberOfPeriods))
    allocate(longitude(size(memberNumber),numberOfPeriods))
    allocate(pressure(size(memberNumber),numberOfPeriods))
    allocate(latitudeWind(size(memberNumber),numberOfPeriods))
    allocate(longitudeWind(size(memberNumber),numberOfPeriods))
    allocate(wind(size(memberNumber),numberOfPeriods))
    allocate(values(size(memberNumber)))
    allocate(period(numberOfPeriods))
    period(1)=0

    ! Observed Storm Centre
    call codes_get(ibufr,'#1#meteorologicalAttributeSignificance',significance);
    call codes_get(ibufr,'#1#latitude',latitudeCentre);
    call codes_get(ibufr,'#1#longitude',longitudeCentre);
    if (significance/=1) then
      print *,'ERROR: unexpected #1#meteorologicalAttributeSignificance'
      stop 1
    endif
    if (latitudeCentre==CODES_MISSING_DOUBLE .and. longitudeCentre==CODES_MISSING_DOUBLE) then
      write(*,'(a)')'Observed storm centre position missing'
    else
      write(*,'(A,F8.2,A,F8.2)')'Observed storm centre: latitude=',latitudeCentre,' longitude=',longitudeCentre
    endif

    ! Location of storm in perturbed analysis
    call codes_get(ibufr,'#2#meteorologicalAttributeSignificance',significance);
    call codes_get(ibufr,'#2#latitude',latitudeAnalysis);
    call codes_get(ibufr,'#2#longitude',longitudeAnalysis);
    call codes_get(ibufr,'#1#pressureReducedToMeanSeaLevel',pressureAnalysis);
    if (significance/=4) then
      print *,'ERROR: unexpected #2#meteorologicalAttributeSignificance'
      stop 1
    endif
    if (size(latitudeAnalysis)==size(memberNumber)) then
      latitude(:,1)=latitudeAnalysis
      longitude(:,1)=longitudeAnalysis
      pressure(:,1)=pressureAnalysis
    else
      latitude(:,1)=latitudeAnalysis(1)
      longitude(:,1)=longitudeAnalysis(1)
      pressure(:,1)=pressureAnalysis(1)
    endif

    ! Location of Maximum Wind
    call codes_get(ibufr,'#3#meteorologicalAttributeSignificance',significance);
    call codes_get(ibufr,'#3#latitude',latitudeMaxWind0);
    call codes_get(ibufr,'#3#longitude',longitudeMaxWind0);
    if (significance/=3) then
      print *,'ERROR: unexpected #3#meteorologicalAttributeSignificance=',significance
      stop 1
    endif
    call codes_get(ibufr,'#1#windSpeedAt10M',windMaxWind0);
    if (size(latitudeMaxWind0)==size(memberNumber)) then
      latitudeWind(:,1)=latitudeMaxWind0
      longitudeWind(:,1)=longitudeMaxWind0
      wind(:,1)=windMaxWind0
    else
      latitudeWind(:,1)=latitudeMaxWind0(1)
      longitudeWind(:,1)=longitudeMaxWind0(1)
      wind(:,1)=windMaxWind0(1)
    endif

    rankSignificance=3
    rankPosition=3
    rankPressure=1
    rankWind=1
    rankPeriod=0

    ! Loop on all periods excluding analysis period(1)=0
    do i=2,numberOfPeriods

      rankPeriod=rankPeriod+1
      write (rankPeriodStr,'(I0)')rankPeriod
      call codes_get(ibufr,'#'//trim(rankPeriodStr)//'#timePeriod',ivalues);
      do k=1,size(ivalues)
        if (ivalues(k)/=CODES_MISSING_LONG) then
          period(i)=ivalues(k)
          exit
        endif
      enddo
      deallocate(ivalues)

      ! Location of the storm
      rankSignificance=rankSignificance+1
      write (rankSignificanceStr,'(I0)')rankSignificance
      call codes_get(ibufr,'#'//trim(rankSignificanceStr)//'#meteorologicalAttributeSignificance',ivalues);
      do k=1,size(ivalues)
        if (ivalues(k)/=CODES_MISSING_LONG) then
          significance=ivalues(k)
          exit
        endif
      enddo
      deallocate(ivalues)

      rankPosition=rankPosition+1
      write (rankPositionStr,'(I0)')rankPosition
      call codes_get(ibufr,'#'//trim(rankPositionStr)//'#latitude',values);
      latitude(:,i)=values
      call codes_get(ibufr,'#'//trim(rankPositionStr)//'#longitude',values);
      longitude(:,i)=values

      if (significance==1) then
        rankPressure=rankPressure+1
        write (rankPressureStr,'(I0)')rankPressure
        call codes_get(ibufr,'#'//trim(rankPressureStr)//'#pressureReducedToMeanSeaLevel',values);
        pressure(:,i)=values
      else
        print *,'ERROR: unexpected meteorologicalAttributeSignificance=',significance
        stop 1
      endif

      ! Location of maximum wind
      rankSignificance=rankSignificance+1
      write (rankSignificanceStr,'(I0)')rankSignificance
      call codes_get(ibufr,'#'//trim(rankSignificanceStr)//'#meteorologicalAttributeSignificance',ivalues);
      do k=1,size(ivalues)
        if (ivalues(k)/=CODES_MISSING_LONG) then
          significance=ivalues(k)
          exit
        endif
      enddo
      deallocate(ivalues)

      rankPosition=rankPosition+1
      write (rankPositionStr,'(I0)')rankPosition
      call codes_get(ibufr,'#'//trim(rankPositionStr)//'#latitude',values);
      latitudeWind(:,i)=values
      call codes_get(ibufr,'#'//trim(rankPositionStr)//'#longitude',values);
      longitudeWind(:,i)=values

      if (significance==3) then
        rankWind=rankWind+1
        write (rankWindStr,'(I0)')rankWind
        call codes_get(ibufr,'#'//trim(rankWindStr)//'#windSpeedAt10M',values);
        wind(:,i)=values
      else
        print *,'ERROR: unexpected meteorologicalAttributeSignificance=,',significance
        stop 1
      endif

    enddo

    ! Print the values
    do i=1,size(memberNumber)
      skipMember=1
      do j=1,size(period)
        if (latitude(i,j)/=CODES_MISSING_DOUBLE .OR. latitudeWind(i,j)/=CODES_MISSING_DOUBLE) then
          skipMember=0
          exit
        endif
      enddo
      if (skipMember/=1) then

        write(*,'(A,I3)') '== Member ',memberNumber(i)
        write(*,*) 'step  latitude   longitude   pressure    latitude   longitude   wind'
        do j=1,size(period)
          if (latitude(i,j)/=CODES_MISSING_DOUBLE .OR. latitudeWind(i,j)/=CODES_MISSING_DOUBLE) then
          write(*,'( I4,2X,F8.2,4X,F8.2,3X,F9.1,2X,F8.2,4X,F8.2,2X,F8.2)') period(j),latitude(i,j),longitude(i,j),pressure(i,j),&
          &latitudeWind(i,j),longitudeWind(i,j),wind(i,j)
          endif
        enddo

      endif
    enddo

    ! deallocating the arrays is very important
    ! because the behaviour of the codes_get functions is as follows:
    ! if the array is not allocated then allocate
    ! if the array is already allocated only copy the values
    deallocate(values)
    deallocate(latitude)
    deallocate(longitude)
    deallocate(pressure)
    deallocate(latitudeWind)
    deallocate(longitudeWind)
    deallocate(wind)
    deallocate(period)
    deallocate(latitudeAnalysis)
    deallocate(longitudeAnalysis)
    deallocate(pressureAnalysis)
    deallocate(memberNumber)
    deallocate(latitudeMaxWind0)
    deallocate(longitudeMaxWind0)
    deallocate(windMaxWind0)

    ! 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_tropical_cyclone

# Copyright 2005-2018 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_tropical_cyclone
#
# Description: how to read data of the ECMWF EPS tropical cyclone tracks encoded in BUFR format.
#
# Please note that tropical cyclone tracks can be encoded in various ways in BUFR.
# Therefore the code below might not work directly for other types of messages
# than the one used in the example. It is advised to use bufr_dump to
# understand the structure of the messages.
#

from __future__ import print_function
import traceback
import sys
import collections

from eccodes import *

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

data = collections.defaultdict(dict)


def example():
    # open BUFR file
    f = open(INPUT, 'rb')

    cnt = 0

    # loop for 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: ', cnt + 1, '  *****************')

        # we need to instruct ecCodes to expand all the descriptors
        # i.e. unpack the data values
        codes_set(bufr, 'unpack', 1)

        numObs = codes_get(bufr, "numberOfSubsets")
        year = codes_get(bufr, "year")
        month = codes_get(bufr, "month")
        day = codes_get(bufr, "day")
        hour = codes_get(bufr, "hour")
        minute = codes_get(bufr, "minute")

        print('Date and time: ', day, '.', month, '.', year, '  ', hour, ':', minute)

        stormIdentifier = codes_get(bufr, "stormIdentifier")
        print('Storm identifier: ', stormIdentifier)

        # How many different timePeriod in the data structure?
        numberOfPeriods = 0
        while True:
            numberOfPeriods = numberOfPeriods + 1
            try:
                codes_get_array(bufr, "#%d#timePeriod" % numberOfPeriods)
            except CodesInternalError as err:
                break
                # the numberOfPeriods includes the analysis (period=0)

        # Get ensembleMemberNumber
        memberNumber = codes_get_array(bufr, "ensembleMemberNumber")
        memberNumberLen = len(memberNumber)

        # Observed Storm Centre
        significance = codes_get(bufr, '#1#meteorologicalAttributeSignificance')
        latitudeCentre = codes_get(bufr, '#1#latitude')
        longitudeCentre = codes_get(bufr, '#1#longitude')

        if significance != 1:
            print('ERROR: unexpected #1#meteorologicalAttributeSignificance')
            return 1

        if (latitudeCentre == CODES_MISSING_DOUBLE) and (longitudeCentre == CODES_MISSING_DOUBLE):
            print('Observed storm centre position missing')
        else:
            print('Observed storm centre: latitude=', latitudeCentre, ' longitude=', longitudeCentre)

        # Location of storm in perturbed analysis
        significance = codes_get(bufr, '#2#meteorologicalAttributeSignificance')

        if significance != 4:
            print('ERROR: unexpected #2#meteorologicalAttributeSignificance')
            return 1

        latitudeAnalysis = codes_get_array(bufr, '#2#latitude')
        longitudeAnalysis = codes_get_array(bufr, '#2#longitude')
        pressureAnalysis = codes_get_array(bufr, '#1#pressureReducedToMeanSeaLevel')

        # Location of Maximum Wind
        significance = codes_get(bufr, '#3#meteorologicalAttributeSignificance')

        if significance != 3:
            print('ERROR: unexpected #3#meteorologicalAttributeSignificance=', significance)
            return 1

        latitudeMaxWind0 = codes_get_array(bufr, '#3#latitude')
        longitudeMaxWind0 = codes_get_array(bufr, '#3#longitude')
        windMaxWind0 = codes_get_array(bufr, '#1#windSpeedAt10M')

        if len(latitudeAnalysis) == len(memberNumber) and len(latitudeMaxWind0) == len(memberNumber):
            for k in range(len(memberNumber)):
                data[k][0] = [latitudeAnalysis[k], longitudeAnalysis[k], pressureAnalysis[k], latitudeMaxWind0[k],
                              longitudeMaxWind0[k], windMaxWind0[k]]

        else:
            for k in range(len(memberNumber)):
                data[k][0] = [latitudeAnalysis[0], longitudeAnalysis[0], pressureAnalysis[k], latitudeMaxWind0[0],
                              longitudeMaxWind0[0], windMaxWind0[k]]

        timePeriod = [0 for x in range(numberOfPeriods)]
        for i in range(1, numberOfPeriods):
            rank1 = i * 2 + 2
            rank3 = i * 2 + 3

            ivalues = codes_get_array(bufr, "#%d#timePeriod" % i)

            if len(ivalues) == 1:
                timePeriod[i] = ivalues[0]
            else:
                for j in range(len(ivalues)):
                    if ivalues[j] != CODES_MISSING_LONG:
                        timePeriod[i] = ivalues[j]
                        break

            # Location of the storm
            values = codes_get_array(bufr, "#%d#meteorologicalAttributeSignificance" % rank1)
            if len(values) == 1:
                significance = values[0]
            else:
                for j in range(len(values)):
                    if values[j] != CODES_MISSING_LONG:
                        significance = values[j]
                        break

            if significance == 1:
                lat = codes_get_array(bufr, "#%d#latitude" % rank1)
                lon = codes_get_array(bufr, "#%d#longitude" % rank1)
                press = codes_get_array(bufr, "#%d#pressureReducedToMeanSeaLevel" % (i + 1))
            else:
                print('ERROR: unexpected meteorologicalAttributeSignificance=', significance)

            # Location of maximum wind
            values = codes_get_array(bufr, "#%d#meteorologicalAttributeSignificance" % rank3)
            if len(values) == 1:
                significanceWind = values[0]
            else:
                for j in range(len(values)):
                    if values[j] != CODES_MISSING_LONG:
                        significanceWind = values[j]
                        break

            if significanceWind == 3:
                latWind = codes_get_array(bufr, "#%d#latitude" % rank3)
                lonWind = codes_get_array(bufr, "#%d#longitude" % rank3)
                wind10m = codes_get_array(bufr, "#%d#windSpeedAt10M" % (i + 1))
            else:
                print('ERROR: unexpected meteorologicalAttributeSignificance=', significanceWind)

            for k in range(len(memberNumber)):
                data[k][i] = [lat[k], lon[k], press[k], latWind[k], lonWind[k], wind10m[k]]


            # ---------------------------------------- Print the values -------------

        for m in range(len(memberNumber)):
            print("== Member  %d" % memberNumber[m])
            print("step  latitude  longitude   pressure  latitude   longitude    wind")
            for s in range(len(timePeriod)):
                if data[m][s][0] != CODES_MISSING_DOUBLE and data[m][s][1] != CODES_MISSING_DOUBLE:
                    print(" {0:>3d}{01}{02:>6.1f}{03}{04:>6.1f}{05}{06:>8.1f}{07}{08:>6.1f}{09}{010:>6.1f}{011}{012:>6.1f}".format(
                        timePeriod[s], '  ', data[m][s][0], '     ', data[m][s][1], '     ', data[m][s][2], '  ',
                        data[m][s][3], '     ', data[m][s][4], '     ', data[m][s][5]))

                # -----------------------------------------------------------------------
        cnt += 1

        # release the BUFR message
        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())