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 as 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 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())