Description
This example shows: how to read data for a tropical cyclone BUFR message.
Source code
bufr_read_tropical_cyclone.f90
! (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.
!
!
! FORTRAN 90 Implementation: bufr_read_tropical_cyclone
!
! Description: how to read data for a tropical cyclone BUFR message.
!
! 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.
!
program bufr_read_tropical_cyclone
use eccodes
implicit none
integer :: ifile
integer :: iret
integer :: ibufr,skipMember
integer :: significance
integer :: year,month,day,hour,minute
integer :: i,j,k,ierr,count=1
integer :: rankPosition,rankSignificance,rankPressure,rankWind
integer :: rankPeriod,numberOfPeriods
real(kind=8) :: latitudeCentre,longitudeCentre
real(kind=8), dimension(:), allocatable :: latitudeMaxWind0,longitudeMaxWind0,windMaxWind0
real(kind=8), dimension(:), allocatable :: latitudeAnalysis,longitudeAnalysis,pressureAnalysis
real(kind=8), dimension(:,:), allocatable :: latitude,longitude,pressure
real(kind=8), dimension(:,:), allocatable :: latitudeWind,longitudeWind,wind
integer(kind=4), dimension(:), allocatable :: memberNumber,period
real(kind=8), dimension(:), allocatable :: values
integer(kind=4), dimension(:), allocatable :: ivalues
character(len=8) :: rankSignificanceStr,rankPositionStr,rankPressureStr,rankWindStr
character(len=8) :: stormIdentifier,rankPeriodStr
call codes_open_file(ifile,'../../data/bufr/tropical_cyclone.bufr','r')
! The first BUFR message is loaded from file.
! ibufr is the BUFR id to be used in subsequent calls
call codes_bufr_new_from_file(ifile,ibufr,iret)
do while (iret/=CODES_END_OF_FILE)
write(*,'(A,I3,A)') '**************** MESSAGE: ',count,' *****************'
! We need to instruct ecCodes to unpack the data values
call codes_set(ibufr,"unpack",1);
call codes_get(ibufr,'year',year);
call codes_get(ibufr,'month',month);
call codes_get(ibufr,'day',day);
call codes_get(ibufr,'hour',hour);
call codes_get(ibufr,'minute',minute);
write(*,'(A,I0,A,I0,A,I0,A,I0,A,I0,A,I0)')'Date and time: ',day,'.',month,'.',year,' ',hour,':',minute
call codes_get(ibufr,'stormIdentifier',stormIdentifier)
write(*,'(A,A)')'Storm identifier: ',stormIdentifier
! How many different timePeriod in the data structure?
rankPeriod=0
ierr=0
do while(ierr==0)
rankPeriod=rankPeriod+1
write (rankPeriodStr,'(I0)')rankPeriod
call codes_get(ibufr,'#'//trim(rankPeriodStr)//'#timePeriod',period,ierr)
if(allocated(period)) deallocate(period)
enddo
! The numberOfPeriods includes the analysis (period=0)
numberOfPeriods=rankPeriod
call codes_get(ibufr,'ensembleMemberNumber',memberNumber)
allocate(latitude(size(memberNumber),numberOfPeriods))
allocate(longitude(size(memberNumber),numberOfPeriods))
allocate(pressure(size(memberNumber),numberOfPeriods))
allocate(latitudeWind(size(memberNumber),numberOfPeriods))
allocate(longitudeWind(size(memberNumber),numberOfPeriods))
allocate(wind(size(memberNumber),numberOfPeriods))
allocate(values(size(memberNumber)))
allocate(period(numberOfPeriods))
period(1)=0
! Observed Storm Centre
call codes_get(ibufr,'#1#meteorologicalAttributeSignificance',significance);
call codes_get(ibufr,'#1#latitude',latitudeCentre);
call codes_get(ibufr,'#1#longitude',longitudeCentre);
if (significance/=1) then
print *,'ERROR: unexpected #1#meteorologicalAttributeSignificance'
stop 1
endif
if (latitudeCentre==CODES_MISSING_DOUBLE .and. longitudeCentre==CODES_MISSING_DOUBLE) then
write(*,'(a)')'Observed storm centre position missing'
else
write(*,'(A,F8.2,A,F8.2)')'Observed storm centre: latitude=',latitudeCentre,' longitude=',longitudeCentre
endif
! Location of storm in perturbed analysis
call codes_get(ibufr,'#2#meteorologicalAttributeSignificance',significance);
call codes_get(ibufr,'#2#latitude',latitudeAnalysis);
call codes_get(ibufr,'#2#longitude',longitudeAnalysis);
call codes_get(ibufr,'#1#pressureReducedToMeanSeaLevel',pressureAnalysis);
if (significance/=4) then
print *,'ERROR: unexpected #2#meteorologicalAttributeSignificance'
stop 1
endif
if (size(latitudeAnalysis)==size(memberNumber)) then
latitude(:,1)=latitudeAnalysis
longitude(:,1)=longitudeAnalysis
pressure(:,1)=pressureAnalysis
else
latitude(:,1)=latitudeAnalysis(1)
longitude(:,1)=longitudeAnalysis(1)
pressure(:,1)=pressureAnalysis(1)
endif
! Location of Maximum Wind
call codes_get(ibufr,'#3#meteorologicalAttributeSignificance',significance);
call codes_get(ibufr,'#3#latitude',latitudeMaxWind0);
call codes_get(ibufr,'#3#longitude',longitudeMaxWind0);
if (significance/=3) then
print *,'ERROR: unexpected #3#meteorologicalAttributeSignificance=',significance
stop 1
endif
call codes_get(ibufr,'#1#windSpeedAt10M',windMaxWind0);
if (size(latitudeMaxWind0)==size(memberNumber)) then
latitudeWind(:,1)=latitudeMaxWind0
longitudeWind(:,1)=longitudeMaxWind0
wind(:,1)=windMaxWind0
else
latitudeWind(:,1)=latitudeMaxWind0(1)
longitudeWind(:,1)=longitudeMaxWind0(1)
wind(:,1)=windMaxWind0(1)
endif
rankSignificance=3
rankPosition=3
rankPressure=1
rankWind=1
rankPeriod=0
! Loop on all periods excluding analysis period(1)=0
do i=2,numberOfPeriods
rankPeriod=rankPeriod+1
write (rankPeriodStr,'(I0)')rankPeriod
call codes_get(ibufr,'#'//trim(rankPeriodStr)//'#timePeriod',ivalues);
do k=1,size(ivalues)
if (ivalues(k)/=CODES_MISSING_LONG) then
period(i)=ivalues(k)
exit
endif
enddo
deallocate(ivalues)
! Location of the storm
rankSignificance=rankSignificance+1
write (rankSignificanceStr,'(I0)')rankSignificance
call codes_get(ibufr,'#'//trim(rankSignificanceStr)//'#meteorologicalAttributeSignificance',ivalues);
do k=1,size(ivalues)
if (ivalues(k)/=CODES_MISSING_LONG) then
significance=ivalues(k)
exit
endif
enddo
deallocate(ivalues)
rankPosition=rankPosition+1
write (rankPositionStr,'(I0)')rankPosition
call codes_get(ibufr,'#'//trim(rankPositionStr)//'#latitude',values);
latitude(:,i)=values
call codes_get(ibufr,'#'//trim(rankPositionStr)//'#longitude',values);
longitude(:,i)=values
if (significance==1) then
rankPressure=rankPressure+1
write (rankPressureStr,'(I0)')rankPressure
call codes_get(ibufr,'#'//trim(rankPressureStr)//'#pressureReducedToMeanSeaLevel',values);
pressure(:,i)=values
else
print *,'ERROR: unexpected meteorologicalAttributeSignificance=',significance
stop 1
endif
! Location of maximum wind
rankSignificance=rankSignificance+1
write (rankSignificanceStr,'(I0)')rankSignificance
call codes_get(ibufr,'#'//trim(rankSignificanceStr)//'#meteorologicalAttributeSignificance',ivalues);
do k=1,size(ivalues)
if (ivalues(k)/=CODES_MISSING_LONG) then
significance=ivalues(k)
exit
endif
enddo
deallocate(ivalues)
rankPosition=rankPosition+1
write (rankPositionStr,'(I0)')rankPosition
call codes_get(ibufr,'#'//trim(rankPositionStr)//'#latitude',values);
latitudeWind(:,i)=values
call codes_get(ibufr,'#'//trim(rankPositionStr)//'#longitude',values);
longitudeWind(:,i)=values
if (significance==3) then
rankWind=rankWind+1
write (rankWindStr,'(I0)')rankWind
call codes_get(ibufr,'#'//trim(rankWindStr)//'#windSpeedAt10M',values);
wind(:,i)=values
else
print *,'ERROR: unexpected meteorologicalAttributeSignificance=,',significance
stop 1
endif
enddo
! Print the values
do i=1,size(memberNumber)
skipMember=1
do j=1,size(period)
if (latitude(i,j)/=CODES_MISSING_DOUBLE .OR. latitudeWind(i,j)/=CODES_MISSING_DOUBLE) then
skipMember=0
exit
endif
enddo
if (skipMember/=1) then
write(*,'(A,I3)') '== Member ',memberNumber(i)
write(*,*) 'step latitude longitude pressure latitude longitude wind'
do j=1,size(period)
if (latitude(i,j)/=CODES_MISSING_DOUBLE .OR. latitudeWind(i,j)/=CODES_MISSING_DOUBLE) then
write(*,'( I4,2X,F8.2,4X,F8.2,3X,F9.1,2X,F8.2,4X,F8.2,2X,F8.2)') period(j),latitude(i,j),longitude(i,j),pressure(i,j),&
&latitudeWind(i,j),longitudeWind(i,j),wind(i,j)
endif
enddo
endif
enddo
! deallocating the arrays is very important
! because the behaviour of the codes_get functions is as follows:
! if the array is not allocated then allocate
! if the array is already allocated only copy the values
deallocate(values)
deallocate(latitude)
deallocate(longitude)
deallocate(pressure)
deallocate(latitudeWind)
deallocate(longitudeWind)
deallocate(wind)
deallocate(period)
deallocate(latitudeAnalysis)
deallocate(longitudeAnalysis)
deallocate(pressureAnalysis)
deallocate(memberNumber)
deallocate(latitudeMaxWind0)
deallocate(longitudeMaxWind0)
deallocate(windMaxWind0)
! Release the BUFR message
call codes_release(ibufr)
! Load the next BUFR message
call codes_bufr_new_from_file(ifile,ibufr,iret)
count=count+1
end do
! Close file
call codes_close_file(ifile)
end program bufr_read_tropical_cyclone
bufr_read_tropical_cyclone.py
# (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())