core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

ECC-196: Update the bufr_read_temp example

This commit is contained in:
Shahram Najm 2015-12-11 14:14:41 +00:00
parent 8bbd1d29fb
commit cfef150e19
5 changed files with 89 additions and 164 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
data/bufr/bufr_data_files.txt
View File

@ -136,3 +136,4 @@ wavb_134.bufr
new.bufr new.bufr
tropical_cyclone.bufr tropical_cyclone.bufr
uegabe.bufr uegabe.bufr
PraticaTemp.bufr

132
examples/F90/bufr_read_temp.f90
View File

@ -1,103 +1,81 @@
! ! Copyright 2005-2015 ECMWF.
!Copyright 2005-2015 ECMWF.
! !
! This software is licensed under the terms of the Apache Licence Version 2.0 ! 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. ! 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 ! 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. ! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
! !
! !
! FOTRAN 90 Implementation: bufr_read_temp ! FORTRAN 90 Implementation: bufr_read_temp
! !
! Description: how to read temperature significant levels from TEMP BUFR messages. ! Description: how to read levels from TEMP BUFR messages.
! !
! Please note that TEMP reports can be encoded in various ways in BUFR. Therefore the code ! 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 ! 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 first to understand the structure of these messages. ! example. It is advised to use bufr_dump first to understand the structure of these messages.
!
program bufr_read_temp program bufr_read_temp
use eccodes use eccodes
implicit none implicit none
integer :: ifile integer :: ifile
integer :: iret integer :: iret,ierr
integer :: ibufr integer :: ibufr
integer :: i, count=0 integer :: i, count=0
integer(kind=4) :: num integer(kind=4),dimension(:), allocatable :: timePeriod,extendedVerticalSoundingSignificance
real(kind=8), dimension(:), allocatable :: presVal,geoVal,tVal ,tdVal integer(kind=4) :: blockNumber,stationNumber,numberOfLevels
character(len=128) :: keyName character(len=30) :: str
real(kind=8),dimension(:), allocatable :: pressure,airTemperature,dewpointTemperature
real(kind=8),dimension(:), allocatable :: geopotentialHeight,latitudeDisplacement,longitudeDisplacement
real(kind=8),dimension(:), allocatable :: windDirection,windSpeed
character(len=128) :: keyName
call codes_open_file(ifile,'../../data/bufr/temp_101.bufr','r') call codes_open_file(ifile,'../../data/bufr/PraticaTemp.bufr','r')
! the first bufr message is loaded from file ! the first bufr message is loaded from file
! ibufr is the bufr id to be used in subsequent calls ! ibufr is the bufr id to be used in subsequent calls
call codes_bufr_new_from_file(ifile,ibufr,iret) call codes_bufr_new_from_file(ifile,ibufr,iret)
do while (iret/=CODES_END_OF_FILE) do while (iret/=CODES_END_OF_FILE)
write(*,*) 'message: ',count write(*,*) 'message: ',count
call codes_set(ibufr,'unpack',1)
! we need to instruct ecCodes to expand all the descriptors call codes_get(ibufr,'timePeriod',timePeriod)
! i.e. unpack the data values call codes_get(ibufr,'pressure',pressure)
call codes_set(ibufr,"unpack",1); call codes_get(ibufr,'extendedVerticalSoundingSignificance',extendedVerticalSoundingSignificance)
call codes_get(ibufr,'geopotentialHeight',geopotentialHeight)
! In what follows we rely on the fact that for call codes_get(ibufr,'latitudeDisplacement',latitudeDisplacement)
! temperature significant levels the value of key call codes_get(ibufr,'longitudeDisplacement',longitudeDisplacement)
! verticalSoundingSignificance is 4 (see flag table 8001 for details). call codes_get(ibufr,'airTemperature',airTemperature)
! call codes_get(ibufr,'dewpointTemperature',dewpointTemperature)
! In our BUFR message verticalSoundingSignificance is always followed by call codes_get(ibufr,'windDirection',windDirection)
! geopotential, airTemperature, dewpointTemperature, call codes_get(ibufr,'windSpeed',windSpeed)
! windDirection, windSpeed and pressure. call codes_get(ibufr,'blockNumber',blockNumber)
! call codes_get(ibufr,'stationNumber',stationNumber)
! So in order to access any of these keys we need to use the print *,'station',blockNumber,stationNumber
! condition: verticalSoundingSignificance=4. print *,'timePeriod pressure geopotentialHeight latitudeDisplacement &
&longitudeDisplacement airTemperature windDirection windSpeed significance'
! ---- Get pressure --------------------------- do i=1,size(windSpeed)
call codes_get(ibufr,'/verticalSoundingSignificance=4/pressure',presVal); write(*,'(I5,6X,F9.1,2X,F9.2,10X,F8.2,14X,F8.2,16X,F8.2,6X,F8.2,4X,F8.2,4X,I0)') timePeriod(i),pressure(i),&
&geopotentialHeight(i),latitudeDisplacement(i),&
! ---- Get gepotential ------------------------ &longitudeDisplacement(i),airTemperature(i),windDirection(i),windSpeed(i),extendedVerticalSoundingSignificance(i)
call codes_get(ibufr,'/verticalSoundingSignificance=4/geopotential',geoVal) enddo
! free arrays
! ---- Get temperature -------------------------------- deallocate(timePeriod)
call codes_get(ibufr,'/verticalSoundingSignificance=4/airTemperature',tVal) deallocate(pressure)
deallocate(geopotentialHeight)
! ---- Get dew point temperature ----------------------- deallocate(latitudeDisplacement)
call codes_get(ibufr,'/verticalSoundingSignificance=4/dewpointTemperature',tdVal) deallocate(longitudeDisplacement)
deallocate(airTemperature)
! ---- Check that all arrays are same size deallocate(dewpointTemperature)
if (size(presVal)/=size(geoVal) .or. size(tVal)/=size(tdVal) .or. size(tdVal)/=size(presVal)) then deallocate(windDirection)
print *,'inconsistent array dimension' deallocate(windSpeed)
exit deallocate(extendedVerticalSoundingSignificance)
endif
num=size(presVal)
! ---- Print the values --------------------------------
write(*,*) 'level pres geo t td'
write(*,*) "--------------------------------------"
do i=1,num
write(*,*) i,presVal(i),geoVal(i),tVal(i),tdVal(i)
end do
! free arrays
deallocate(presVal)
deallocate(geoVal)
deallocate(tVal)
deallocate(tdVal)
! release the bufr message ! release the bufr message
call codes_release(ibufr) call codes_release(ibufr)
! load the next bufr message ! load the next bufr message
call codes_bufr_new_from_file(ifile,ibufr,iret) call codes_bufr_new_from_file(ifile,ibufr,iret)
count=count+1 count=count+1
end do
end do ! close file
! close file
call codes_close_file(ifile) call codes_close_file(ifile)
end program bufr_read_temp end program bufr_read_temp

9
examples/F90/bufr_read_temp.sh
View File

@ -18,18 +18,15 @@ label="bufr_read_temp_f"
fTmp=${label}.tmp.txt fTmp=${label}.tmp.txt
rm -f $fTmp | true rm -f $fTmp | true
#We check "temp_101.bufr". The path is # The path to the BUFR file is hardcoded in the example
#hardcoded in the example
REDIRECT=/dev/null REDIRECT=/dev/null
#Write the key values into a file # Run the example
${examples_dir}/eccodes_f_bufr_read_temp 2> $REDIRECT > $fTmp ${examples_dir}/eccodes_f_bufr_read_temp 2> $REDIRECT > $fTmp
#TODO: check the results #TODO: check the results
#cat $fTmp
#Clean up #Clean up
rm -f $fTmp | true rm -f $fTmp

101
examples/python/bufr_read_temp.py
View File

@ -1,3 +1,4 @@
#
# Copyright 2005-2015 ECMWF. # Copyright 2005-2015 ECMWF.
# #
# This software is licensed under the terms of the Apache Licence Version 2.0 # This software is licensed under the terms of the Apache Licence Version 2.0
@ -6,12 +7,11 @@
# In applying this licence, ECMWF does not waive the privileges and immunities # In applying this licence, ECMWF does not waive the privileges and immunities
# granted to it by virtue of its status as an intergovernmental organisation # granted to it by virtue of its status as an intergovernmental organisation
# nor does it submit to any jurisdiction. # nor does it submit to any jurisdiction.
# #
# Python implementation: bufr_read_temp # Python implementation: bufr_read_temp
# #
# Description: how to read temperature significant levels from TEMP BUFR # Description: how to read levels from TEMP BUFR messages.
# messages. #
# #
# #
# Please note that TEMP reports can be encoded in various ways in BUFR. # Please note that TEMP reports can be encoded in various ways in BUFR.
@ -19,94 +19,50 @@
# messages than the one used in the example. It is advised to use bufr_dump to # messages than the one used in the example. It is advised to use bufr_dump to
# understand the structure of the messages. # understand the structure of the messages.
# #
import traceback import traceback
import sys import sys
from eccodes import * from eccodes import *
INPUT = '../../data/bufr/temp_101.bufr' INPUT = '../../data/bufr/PraticaTemp.bufr'
VERBOSE = 1 # verbose error reporting VERBOSE = 1 # verbose error reporting
def example(): def example():
# open bufr file # open bufr file
f = open(INPUT) f = open(INPUT)
cnt = 0 cnt = 0
# loop over the messages in the file
# loop for the messages in the file
while 1: while 1:
# get handle for message # get handle for message
gid = codes_bufr_new_from_file(f) bufr = codes_bufr_new_from_file(f)
if gid is None: if bufr is None:
break break
print "message: %s" % cnt print "message: %s" % cnt
# we need to instruct ecCodes to expand all the descriptors # we need to instruct ecCodes to expand all the descriptors
# i.e. unpack the data values # i.e. unpack the data section
codes_set(gid, 'unpack', 1) codes_set(bufr, 'unpack', 1)
# get all the timePeriods
# In what follows we rely on the fact that for timePeriod = codes_get_array(bufr, "timePeriod")
# temperature significant levels the value of key pressure = codes_get_array(bufr, "pressure")
# verticalSoundingSignificance is 4 (see flag table 8001 for details). extendedVerticalSoundingSignificance = codes_get_array(bufr, "extendedVerticalSoundingSignificance")
geopotentialHeight = codes_get_array(bufr, "geopotentialHeight")
# We also make use of the fact that in our BUFR message latitudeDisplacement = codes_get_array(bufr, "latitudeDisplacement")
# verticalSoundingSignificance is always followed by geopotential, longitudeDisplacement = codes_get_array(bufr, "longitudeDisplacement")
# airTemperature, dewpointTemperature, airTemperature = codes_get_array(bufr, "airTemperature")
# windDirection, windSpeed and pressure. dewpointTemperature = codes_get_array(bufr, "dewpointTemperature")
windDirection = codes_get_array(bufr, "windDirection")
# Get the number of the temperature significant levels. windSpeed = codes_get_array(bufr, "windSpeed")
blockNumber = codes_get(bufr, "blockNumber")
# We find out the number of temperature significant levels by stationNumber = codes_get(bufr, "stationNumber")
# counting how many pressure values we have on these levels. print 'station %d%d' % (blockNumber,stationNumber)
print 'timePeriod pressure geopotentialHeight latitudeDisplacement longitudeDisplacement airTemperature windDirection windSpeed significance'
numSigT = codes_get_size( for i in range(0,len(windSpeed)-1):
gid, "/verticalSoundingSignificance=4/pressure") print timePeriod[i],pressure[i],geopotentialHeight[i],latitudeDisplacement[i],longitudeDisplacement[i],airTemperature[i],windDirection[i],windSpeed[i],extendedVerticalSoundingSignificance[i]
print ' Number of temperature significant levels %ld' % (numSigT)
# Get pressure
sigt_pres = codes_get_array(
gid, "/verticalSoundingSignificance=4/pressure")
# Get gepotential
sigt_geo = codes_get_array(
gid, "/verticalSoundingSignificance=4/geopotential")
# Get temperature
sigt_t = codes_get_array(
gid, "/verticalSoundingSignificance=4/airTemperature")
# Get dew point
sigt_td = codes_get_array(
gid, "/verticalSoundingSignificance=4/dewpointTemperature")
# Check that all arrays are same size
if len(sigt_pres) != numSigT or len(sigt_geo) != numSigT or \
len(sigt_t) != numSigT or len(sigt_td) != numSigT:
print 'inconsistent array dimension'
return 1
# Print the values
print "lev pres geo t td"
print "-------------------------------"
for i in xrange(numSigT):
print "%3d %6.0f %6.0f %.1f %.1f" % (i + 1, sigt_pres[i],
sigt_geo[i], sigt_t[i],
sigt_td[i])
cnt += 1 cnt += 1
# delete handle # delete handle
codes_release(gid) codes_release(bufr)
# close the file # close the file
f.close() f.close()
def main(): def main():
try: try:
example() example()
@ -115,8 +71,7 @@ def main():
traceback.print_exc(file=sys.stderr) traceback.print_exc(file=sys.stderr)
else: else:
print >>sys.stderr, err.msg print >>sys.stderr, err.msg
return 1 return 1
if __name__ == "__main__": if __name__ == "__main__":
sys.exit(main()) sys.exit(main())

10
examples/python/bufr_read_temp.sh
View File

@ -16,20 +16,14 @@ label="bufr_read_temp_p"
#Define tmp file #Define tmp file
fTmp=${label}.tmp.txt fTmp=${label}.tmp.txt
rm -f $fTmp | true rm -f $fTmp
#We check "temp_101.bufr". The path is
#hardcoded in the example
REDIRECT=/dev/null REDIRECT=/dev/null
#Write the key values into a file # Run it. The path to the BUFR file is hardcoded in the example
$PYTHON $examples_src/bufr_read_temp.py 2> $REDIRECT > $fTmp $PYTHON $examples_src/bufr_read_temp.py 2> $REDIRECT > $fTmp
#TODO: check the results #TODO: check the results
#cat $fTmp
#Clean up #Clean up
rm -f $fTmp rm -f $fTmp