mirror of https://github.com/ecmwf/eccodes.git
268 lines
9.0 KiB
Python
268 lines
9.0 KiB
Python
# (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())
|