eccodes/examples/python/bufr_read_tropical_cyclone.py

# Copyright 2005-2016 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.
#

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)

    cnt = 0

    # loop for the messages in the file
    while 1:
        # get handle for message
        gid = codes_bufr_new_from_file(f)
        if gid 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(gid, 'unpack', 1)

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

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

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

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

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

        # Observed Storm Centre
        significance    = codes_get(gid,'#1#meteorologicalAttributeSignificance')
        latitudeCentre  = codes_get(gid,'#1#latitude')
        longitudeCentre = codes_get(gid,'#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(gid,'#2#meteorologicalAttributeSignificance')

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

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

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

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

        latitudeMaxWind0=codes_get_array(gid,'#3#latitude')
        longitudeMaxWind0= codes_get_array(gid,'#3#longitude')
        windMaxWind0= codes_get_array(gid,'#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(gid,"#%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(gid, "#%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(gid, "#%d#latitude" % rank1)
                lon = codes_get_array(gid, "#%d#longitude" % rank1)
                press = codes_get_array(gid, "#%d#pressureReducedToMeanSeaLevel" % (i + 1))
            else:
                print 'ERROR: unexpected meteorologicalAttributeSignificance=',significance

            #Location of maximum wind
            values = codes_get_array(gid, "#%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(gid, "#%d#latitude" % rank3)
                lonWind = codes_get_array(gid, "#%d#longitude" % rank3)
                wind10m = codes_get_array(gid, "#%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 " {:>3d}{}{:>6.1f}{}{:>6.1f}{}{:>8.1f}{}{:>6.1f}{}{:>6.1f}{}{:>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(gid)

    # 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())