# (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_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. # import collections import sys import traceback 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}{1}{2:>6.1f}{3}{4:>6.1f}{5}{6:>8.1f}{7}{8:>6.1f}{9}{10:>6.1f}{11}{12:>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())