Description
This example shows: how to read data for a tropical cyclone BUFR message.Source code
bufr_read_tropical_cyclone.f90
!Copyright 2005-2016 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.
!
!
! FOTRAN 90 Implementation: bufr_read_tropical_cyclone
!
! Description: how to read data for a tropical cyclone BUFR message.
!
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
# Copyright 2005-2016 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.
#
import traceback
import sys
import collections
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)
cnt = 0
# loop for the messages in the file
while 1:
# get handle for message
gid = codes_bufr_new_from_file(f)
if gid 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(gid, 'unpack', 1)
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: ', day,'.',month,'.',year,' ',hour,':',minute
stormIdentifier = codes_get(gid,"stormIdentifier")
print 'Storm identifier: ', stormIdentifier
#How many different timePeriod in the data structure?
numberOfPeriods=0
while True:
numberOfPeriods=numberOfPeriods+1
try:
codes_get_array(gid,"#%d#timePeriod" %numberOfPeriods)
except CodesInternalError, err:
break
#the numberOfPeriods includes the analysis (period=0)
# Get ensembleMemberNumber
memberNumber = codes_get_array(gid, "ensembleMemberNumber")
memberNumberLen=len(memberNumber)
# Observed Storm Centre
significance = codes_get(gid,'#1#meteorologicalAttributeSignificance')
latitudeCentre = codes_get(gid,'#1#latitude')
longitudeCentre = codes_get(gid,'#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(gid,'#2#meteorologicalAttributeSignificance')
if significance!=4:
print 'ERROR: unexpected #2#meteorologicalAttributeSignificance'
return 1
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')
if significance!=3:
print 'ERROR: unexpected #3#meteorologicalAttributeSignificance=', significance
return 1
latitudeMaxWind0=codes_get_array(gid,'#3#latitude')
longitudeMaxWind0= codes_get_array(gid,'#3#longitude')
windMaxWind0= codes_get_array(gid,'#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(gid,"#%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(gid, "#%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(gid, "#%d#latitude" % rank1)
lon = codes_get_array(gid, "#%d#longitude" % rank1)
press = codes_get_array(gid, "#%d#pressureReducedToMeanSeaLevel" % (i + 1))
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]
else:
for j in range(len (values)):
if values[j]!=CODES_MISSING_LONG:
significanceWind=values[j]
break
if significanceWind==3:
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(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 " {:>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)
# close the file
f.close()
def main():
try:
example()
except CodesInternalError, err:
if VERBOSE:
traceback.print_exc(file=sys.stderr)
else:
print >>sys.stderr, err.msg
return 1
if __name__ == "__main__":
sys.exit(main())