mirror of https://github.com/ecmwf/eccodes.git
tropical cyclone example in python ECC-190
This commit is contained in:
Marijana Crepulja
parent
b049f81f84
commit
7a541afa92
|
|
@ -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()
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue