mirror of https://github.com/ecmwf/eccodes.git
225 lines
7.7 KiB
Python
225 lines
7.7 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_tempf
|
|
#
|
|
#
|
|
# Description: read and print radiosonde data from TEMP BUFR messages.
|
|
# If available this version also lists the position information from the WMO list
|
|
# (now OSCAR/Surface) appended to the reports by ECMWF
|
|
#
|
|
# Author: Bruce Ingleby
|
|
#
|
|
# Please note that TEMP reports can be encoded in various ways in BUFR.
|
|
# Therefore the code below might not work directly for other types of TEMP
|
|
# messages than the one used in the example. It is advised to use bufr_dump to
|
|
# understand the structure of the messages.
|
|
#
|
|
|
|
import sys
|
|
import traceback
|
|
|
|
import numpy as np
|
|
|
|
from eccodes import *
|
|
|
|
INPUT = "../../data/bufr/PraticaTemp.bufr"
|
|
VERBOSE = 1 # verbose error reporting
|
|
|
|
|
|
def example():
|
|
# open BUFR file
|
|
f = open(INPUT, "rb")
|
|
llstdonly = 1 # If 1 then list standard levels only, 0 list all levels
|
|
cnt = 0
|
|
# loop over the messages in the file
|
|
while 1:
|
|
# get handle for message
|
|
bufr = codes_bufr_new_from_file(f)
|
|
if bufr is None:
|
|
break
|
|
cnt += 1
|
|
|
|
# desc = codes_get_array(bufr, 'unexpandedDescriptors')
|
|
# if all(desc != 309056): # descent reports
|
|
# codes_release(bufr)
|
|
# continue # Skip other templates
|
|
|
|
# we need to instruct ecCodes to expand all the descriptors
|
|
# i.e. unpack the data section
|
|
codes_set(bufr, "unpack", 1)
|
|
# get header information from the message
|
|
try:
|
|
dropid = codes_get(bufr, "aircraftRegistrationNumberOrOtherIdentification")
|
|
except Exception:
|
|
dropid = "UNKNOWN "
|
|
try:
|
|
shipid = codes_get(bufr, "shipOrMobileLandStationIdentifier")
|
|
except Exception:
|
|
shipid = "UNKNOWN "
|
|
|
|
statid = "00000 "
|
|
try:
|
|
block = codes_get(bufr, "blockNumber")
|
|
stnum = codes_get(bufr, "stationNumber")
|
|
if (block > 0) and (block < 100): # or block != CODES_MISSING_LONG
|
|
statid = str.format("%.2i%.3i " % (block, stnum))
|
|
except Exception:
|
|
statid = "00000 "
|
|
if statid == "00000 ":
|
|
statid = shipid[0:8]
|
|
if statid == "UNKNOWN ":
|
|
statid = dropid[0:8]
|
|
|
|
# subtype = codes_get(bufr,'rdbSubtype')
|
|
sondetype = codes_get(bufr, "radiosondeType")
|
|
if sondetype == CODES_MISSING_LONG:
|
|
sondetype = 0
|
|
slat = codes_get_array(bufr, "latitude")
|
|
slon = codes_get_array(bufr, "longitude")
|
|
slat = np.where(slat != CODES_MISSING_DOUBLE, slat, np.nan)
|
|
slon = np.where(slon != CODES_MISSING_DOUBLE, slon, np.nan)
|
|
try:
|
|
htg = codes_get(bufr, "heightOfStationGroundAboveMeanSeaLevel")
|
|
if htg == CODES_MISSING_DOUBLE:
|
|
htg = np.nan
|
|
except Exception:
|
|
htg = np.nan
|
|
try:
|
|
htp = codes_get(bufr, "heightOfBarometerAboveMeanSeaLevel")
|
|
if htp == CODES_MISSING_DOUBLE:
|
|
htp = np.nan
|
|
except Exception:
|
|
htp = np.nan
|
|
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")
|
|
try:
|
|
second = codes_get(bufr, "second")
|
|
if second == CODES_MISSING_LONG:
|
|
second = 0
|
|
except Exception:
|
|
second = 0
|
|
date = str.format("%i%.2i%.2i" % (year, month, day))
|
|
time = str.format("%.2i%.2i%.2i" % (hour, minute, second))
|
|
try:
|
|
windsp = codes_get_array(bufr, "windSpeed")
|
|
except Exception:
|
|
codes_release(bufr)
|
|
continue
|
|
print(
|
|
"Ob: %7i %8s %s %s %7.3f %8.3f %6.1f %6.1f %3i %4i"
|
|
% (
|
|
cnt,
|
|
statid,
|
|
date,
|
|
time,
|
|
slat[0],
|
|
slon[0],
|
|
htg,
|
|
htp,
|
|
sondetype,
|
|
len(windsp),
|
|
)
|
|
)
|
|
|
|
try:
|
|
rsnumber = codes_get(bufr, "radiosondeSerialNumber")
|
|
rssoftware = codes_get(bufr, "softwareVersionNumber")
|
|
balloonwt = codes_get(bufr, "weightOfBalloon")
|
|
print("RS number/software/balloonwt ", rsnumber, rssoftware, balloonwt)
|
|
except Exception:
|
|
rsnumber = 0
|
|
try:
|
|
htec = codes_get(
|
|
bufr, "heightOfStation"
|
|
) # Height from WMO list (appended by ECMWF)
|
|
if htec == CODES_MISSING_LONG:
|
|
htec = np.nan
|
|
print(
|
|
"WMO list lat, lon, ht: %s %7.3f %8.3f %6.1f"
|
|
% (statid, slat[1], slon[1], htec)
|
|
)
|
|
except Exception:
|
|
htec = np.nan
|
|
|
|
# get all the timePeriods
|
|
dtime = codes_get_array(bufr, "timePeriod")
|
|
try:
|
|
pressure = codes_get_array(bufr, "pressure")
|
|
except Exception:
|
|
codes_release(bufr)
|
|
continue
|
|
vsSignif = codes_get_array(bufr, "extendedVerticalSoundingSignificance")
|
|
try:
|
|
geopoth = codes_get_array(bufr, "nonCoordinateGeopotentialHeight")
|
|
except Exception:
|
|
codes_release(bufr)
|
|
continue
|
|
dlat = codes_get_array(bufr, "latitudeDisplacement")
|
|
dlon = codes_get_array(bufr, "longitudeDisplacement")
|
|
airt = codes_get_array(bufr, "airTemperature")
|
|
dewt = codes_get_array(bufr, "dewpointTemperature")
|
|
windd = codes_get_array(bufr, "windDirection")
|
|
dtime = np.where(dtime != CODES_MISSING_LONG, dtime, np.nan)
|
|
dlat = np.where(dlat != CODES_MISSING_DOUBLE, dlat, np.nan)
|
|
dlon = np.where(dlon != CODES_MISSING_DOUBLE, dlon, np.nan)
|
|
airt = np.where(airt != CODES_MISSING_DOUBLE, airt, np.nan)
|
|
dewt = np.where(dewt != CODES_MISSING_DOUBLE, dewt, np.nan)
|
|
windd = np.where(windd != CODES_MISSING_LONG, windd, np.nan)
|
|
windsp = np.where(windsp != CODES_MISSING_DOUBLE, windsp, np.nan)
|
|
geopoth = np.where(geopoth != CODES_MISSING_LONG, geopoth, np.nan)
|
|
pressure = np.where(pressure != CODES_MISSING_DOUBLE, pressure, np.nan)
|
|
# pressure = np.where(pressure > -1e10, pressure, np.nan)
|
|
print(
|
|
"level dtime dlat dlon pressure geopotH airTemp dewPtT windDir windSp signif"
|
|
)
|
|
for i in range(0, len(windsp)):
|
|
if llstdonly and vsSignif[i] != 65536:
|
|
continue
|
|
print(
|
|
"%5i %6.1f %6.3f %6.3f %8.1f %7.1f %7.2f %7.2f %7.2f %7.2f %7i"
|
|
% (
|
|
i + 1,
|
|
dtime[i],
|
|
dlat[i],
|
|
dlon[i],
|
|
pressure[i],
|
|
geopoth[i],
|
|
airt[i],
|
|
dewt[i],
|
|
windd[i],
|
|
windsp[i],
|
|
vsSignif[i],
|
|
)
|
|
)
|
|
# delete handle
|
|
codes_release(bufr)
|
|
# close the file
|
|
f.close()
|
|
print("Finishing normally. Number of BUFR records read: ", cnt)
|
|
|
|
|
|
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())
|