From 7a541afa925552681f8c177b3ac7753cd5f4b72b Mon Sep 17 00:00:00 2001 From: Marijana Crepulja Date: Wed, 23 Dec 2015 15:25:57 +0000 Subject: [PATCH] tropical cyclone example in python ECC-190 --- examples/python/bufr_read_tropical_cyclone.py | 150 +++++------------- 1 file changed, 40 insertions(+), 110 deletions(-) diff --git a/examples/python/bufr_read_tropical_cyclone.py b/examples/python/bufr_read_tropical_cyclone.py index 8615dcb01..2881905be 100644 --- a/examples/python/bufr_read_tropical_cyclone.py +++ b/examples/python/bufr_read_tropical_cyclone.py @@ -10,16 +10,10 @@ # # Description: how to read data of the ECMWF EPS tropical cyclone tracks encoded in BUFR format. # -# Please note that scatterometer data 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 first to understand the structure of these messages. -# - import traceback import sys import collections -#import numpy from eccodes import * @@ -28,12 +22,13 @@ VERBOSE = 1 # verbose error reporting data=collections.defaultdict(dict) + def example(): # open bufr file f = open(INPUT) - cnt = 1 + cnt = 0 # loop for the messages in the file while 1: @@ -42,52 +37,37 @@ def example(): if gid is None: break - print '**************** MESSAGE: ',cnt,' *****************' + 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) -# The BUFR file contains a single message with number of subsets in a compressed form. -# It means each subset has exactly the same structure. -# One subset contains values of latitude, longitude, pressure, wind at 10m -# for one particular ensemble member over forecast period from 0h to 240h by 6h step. -# To print values of latitude, longitude, pressure, wind at 10m for particular ensemble member from all the subsets -# we will simply access the key by condition (see below) -# -# The key latitude will give back the array of all the values corresponding -# to all the instances of the key in the tree. -# The key #2#latitude will return only the values of the second instance in the tree. - - numObs=codes_get(gid,"numberOfSubsets") + 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: %d%s %d%s %d%s %d%s %d' % (day,'.',month,'.',year,' ',hour,':',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? - rankPeriod=0 + numberOfPeriods=0 while True: - rankPeriod=rankPeriod+1 + numberOfPeriods=numberOfPeriods+1 try: - codes_get_array(gid,"#%d#timePeriod" %rankPeriod) - #period=next(timeStep[i] for i in range (len(timeStep))if timeStep[i] != CODES_MISSING_DOUBLE ) + codes_get_array(gid,"#%d#timePeriod" %numberOfPeriods) except CodesInternalError, err: - #print err.msg - #del(period) break #the numberOfPeriods includes the analysis (period=0) + # Get ensembleMemberNumber memberNumber = codes_get_array(gid, "ensembleMemberNumber") memberNumberLen=len(memberNumber) @@ -109,72 +89,40 @@ def example(): # Location of storm in perturbed analysis significance = codes_get(gid,'#2#meteorologicalAttributeSignificance') - latitudeAnalysis = codes_get_array(gid,'#2#latitude') - longitudeAnalysis = codes_get_array(gid,'#2#longitude') - pressureAnalysis = codes_get_array(gid,'#1#pressureReducedToMeanSeaLevel') if (significance!=4): print 'ERROR: unexpected #2#meteorologicalAttributeSignificance' return 1 - latitude=[[0 for x in range(rankPeriod)] for x in range(len(memberNumber))] - longitude=[[0 for x in range(rankPeriod)] for x in range(len(memberNumber))] - pressure=[[0 for x in range(rankPeriod)] for x in range(len(memberNumber))] - - - if (len(latitudeAnalysis)==len(memberNumber)): - for i in range(len(memberNumber)): - latitude[i][0]=latitudeAnalysis[i] - longitude[i][0]=longitudeAnalysis[i] - pressure[i][0]=pressureAnalysis[i] - - else: - for i in range(len(memberNumber)): - latitude[i][0]=latitudeAnalysis[0] - longitude[i][0]=longitudeAnalysis[0] - pressure[i][0]=pressureAnalysis[i] + 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') - latitudeMaxWind0=codes_get_array(gid,'#3#latitude') - longitudeMaxWind0= codes_get_array(gid,'#3#longitude') if (significance!=3): print 'ERROR: unexpected #3#meteorologicalAttributeSignificance=', significance return 1 - latitudeWind=[[0 for x in range(rankPeriod)] for x in range(len(memberNumber))] - longitudeWind=[[0 for x in range(rankPeriod)] for x in range(len(memberNumber))] - wind=[[0 for x in range(rankPeriod)] for x in range(len(memberNumber))] - + latitudeMaxWind0=codes_get_array(gid,'#3#latitude') + longitudeMaxWind0= codes_get_array(gid,'#3#longitude') windMaxWind0= codes_get_array(gid,'#1#windSpeedAt10M') - if (len(latitudeMaxWind0)==len(memberNumber)): - for i in range(len(memberNumber)): - latitudeWind[i][0]=latitudeMaxWind0[i] - longitudeWind[i][0]=longitudeMaxWind0[i] - wind[i][0]=windMaxWind0[i] + 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 i in range(len(memberNumber)): - latitudeWind[i][0]=latitudeMaxWind0[0] - longitudeWind[i][0]=longitudeMaxWind0[0] - wind[i][0]=windMaxWind0[i] + for k in range(len(memberNumber)): + data[k][0]=[latitudeAnalysis[0],longitudeAnalysis[0],pressureAnalysis[k],latitudeMaxWind0[0],longitudeMaxWind0[0],windMaxWind0[k]] - rankSignificance=3 - rankPosition=3 - rankPressure=1 - rankWind=1 - timePeriod=[0 for x in range(rankPeriod)] - for i in range(1,rankPeriod): - #for i in range(rankPeriod): + timePeriod=[0 for x in range(numberOfPeriods)] + for i in range(1,numberOfPeriods): rank1 = i * 2 + 2 rank3 = i * 2 + 3 - #rank1 = (i+1) * 2 + 2 - #rank3 = (i+1)* 2 + 3 - ivalues= codes_get_array(gid,"#%d#timePeriod" %(i)) if (len(ivalues)==1): @@ -185,6 +133,7 @@ def example(): timePeriod[i]=ivalues[j] break + #Location of the storm values = codes_get_array(gid, "#%d#meteorologicalAttributeSignificance" % rank1) if (len(values)==1): significance=values[0] @@ -193,7 +142,7 @@ def example(): 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) @@ -201,7 +150,7 @@ def example(): 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] @@ -212,54 +161,35 @@ def example(): break if(significanceWind==3): - lat3 = codes_get_array(gid, "#%d#latitude" % rank3) - lon3 = codes_get_array(gid, "#%d#longitude" % rank3) - wind3 = codes_get_array(gid, "#%d#windSpeedAt10M" % (i + 1)) + 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(1,len(memberNumber)): + for k in range(len(memberNumber)): - data[k][i]=[lat[k],lon[k],press[k],lat3[k],lon3[k],wind3[k]] + data[k][i]=[lat[k],lon[k],press[k],latWind[k],lonWind[k],wind10m[k]] - if(cnt==2): - print 'majaaa', i,timePeriod[s] - # ---- Print the values -------------------------------- +# ---------------------------------------- Print the values ----------------------------------------------- - for m in sorted( data.keys()): - print "== Member %d"%(m+1) + for m in range(len(memberNumber)): + print "== Member %d" %memberNumber[m] print "step latitude longitude pressure latitude longitude wind" - #print " {:<3d}{:2.2f}{:2.2f}{}{} {} ".format(timePeriod[0],latitude[m][0],longitude[m][0],pressure[m][0],latitudeWind[m][0],longitudeWind[i][0],wind[i][0]) - if latitude[m][0]!=CODES_MISSING_DOUBLE: - print timePeriod[0],latitude[m][0],longitude[m][0],pressure[m][0], - #print " {:>3d}{}{:>6.1f}{}{:>6.1f}{}{:>8.1f}{}{:>6.1f}{}{:>6.1f}{}{:>6.1f}".format(\ - # timePeriod[0],' ',latitude[m][0],' ',longitude[m][0],' ',pressure[m][0],' ',latitudeWind[m][0],' ',longitudeWind[m][0],' ',wind[m][0]) - - for s in sorted( data[m].keys()): - if data[m][s][0]!=CODES_MISSING_DOUBLE: - print 'maja',s - #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]) - #print "%d %.6g %.6g"%(timePeriod[s],data[m][s][0],data[m][s][1]) - + 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) - cnt += 1 - print 'majaaaaaaaaaaa',cnt - - #del(timePeriod) - #del(latitude) - #del(longitude) - #del(pressure) - #del(latitudeWind) - #del(longitudeWind) - #del(wind) - #del(data) - - # close the file + # close the file f.close()