使用观测资料编码prepbufr文件（python）

prepbufr编码

prepbufr数据的解码在之前的文章有介绍过，这里简单介绍一下利用已知的观测数据（excel格式）取编码prepbufr二进制数据

---

---

一、编码

ncepbufr的git网页给的编码程序text_write.py
值得注意的是，它只适用于添加或者补录单个观测数据,
我的目的是循环读取不同时间、层次、属性、类型数据并编码进去。

from __future__ import print_function
import ncepbufr
import numpy as np

# run test.py first to create prepbufr.table file.

hdstr='SID XOB YOB DHR TYP ELV SAID T29'
obstr='POB QOB TOB ZOB UOB VOB PWO MXGS HOVI CAT PRSS TDO PMO'
qcstr='PQM QQM TQM ZQM WQM NUL PWQ PMQ'
oestr='POE QOE TOE NUL WOE NUL PWE'

# open prepbufr file.

bufr = ncepbufr.open('data/prepbufr2','w',table='prepbufr.table')
idate=2010050700 # cycle time: YYYYMMDDHH
subset='ADPSFC'  # surface land (SYNOPTIC, METAR) reports
bufr.open_message(subset, idate)

hdr = bufr.missing_value*np.ones(len(hdstr.split()),np.float)
hdr[0] = np.fromstring('KTKI    ',dtype=np.float)[0]
hdr[1]=263.4; hdr[2]=33.2; hdr[3] = -0.1; hdr[4]=287; hdr[5]=179
# encode header for wind obs
bufr.write_subset(hdr,hdstr)
# set obs, qcf, oer for  wind
obs = bufr.missing_value*np.ones(len(obstr.split()),np.float)
oer = bufr.missing_value*np.ones(len(oestr.split()),np.float)
qcf = bufr.missing_value*np.ones(len(qcstr.split()),np.float)
obs[0]=985.2; obs[4]=-2.8; obs[5]=-7.7; obs[7]=6.0
qcf[0]=2.0  ; qcf[4]=2.0; oer[4] = 1.6
# encode wind obs
bufr.write_subset(obs,obstr)
# encode wind ob err
bufr.write_subset(oer,oestr)
# encode quality flags, end subset.
bufr.write_subset(qcf,qcstr,end=True)
# set obs, qcf, oer for  temperature and moisture
hdr[4]=187          # report type
# encode header
bufr.write_subset(hdr,hdstr)
obs[:]=bufr.missing_value; qcf[:]=bufr.missing_value; oer[:]=bufr.missing_value
obs[0]=985.2;obs[1]=12968.0;obs[2]=31.3;obs[3]=179.0;obs[7]=0.0
qcf[0]=2.0  ;qcf[1]=2.0    ;qcf[2]=2.0 ;qcf[3]=2.0
oer[0]=0.5  ;oer[1]=0.6    ;oer[2]=2.3
# encode temperature and moisture obs, obs err and qc flags.
bufr.write_subset(obs,obstr)
bufr.write_subset(oer,oestr)
bufr.write_subset(qcf,qcstr,end=True) # end subset
# close bufr message
bufr.close_message()
# create a message with upper-air (raob) data.
subset='ADPUPA'  # upper-air (raob, drops) reports
bufr.open_message(subset, idate)

# set header
hdr[:]=bufr.missing_value
hdr[0] = np.fromstring('72293   ',dtype=np.float)[0]
hdr[1]=242.9; hdr[2]=32.9; hdr[3]=0.0; hdr[5]=134.0

# set obs, qcf, oer for  wind
nlvl=3
hdr[4]=220          # report type: sounding
obs = bufr.missing_value*np.ones