IODA Bundle
JEDI_amsu_bufr2nc.v8.py
Go to the documentation of this file.
1 #!/usr/bin/env python3
2 """convert amsubufr file to netCDF"""
3 from __future__ import print_function
4 import ncepbufr
5 import numpy as np
6 import netCDF4
7 import numbers
8 import sys
9 from netCDF4 import Dataset
10 import re
11 
12 # original routine from Scott Gregory
13 # CSD: cleaned up, matched to gsi obs_diag output
14 # Todo: test for presence of CLATH, CLONH
15 
16 if len(sys.argv) < 3:
17  raise SystemExit('radbufr2nc <input prepbufr> <output netcdf PATH>')
18 ##radbufr_filename='1bamua'
19 # input and output file names from command line args.
20 radbufr_filename = sys.argv[1]
21 netcdf_filePATH=sys.argv[2]
22 
23 ## EXAMPLE BUFR /pan2/projects/gfsenkf/whitaker/cfsr_dumps/2005022506/gdas/1bamua.gdas.2005022506
24 datefrombufr=re.findall(r'(\d{10})', radbufr_filename) #finds ten digit number in the string preceded by PERIOD... ## EXAMPLE BUFR /pan2/projects/gfsenkf/whitaker/cfsr_dumps/2005022506/gdas/1bamua.gdas.2005022506
25 try:
26  filetime=datefrombufr[len(datefrombufr)-1]
27 except:
28  filetime=datefrombufr
29 
30 ##datefrombufr=re.findall(r'(\d{10})', radbufr_filename) #finds ten digit number in the string preceded by PERIOD... ## EXAMPLE BUFR /pan2/projects/gfsenkf/whitaker/cfsr_dumps/2005022506/gdas/1bamua.gdas.2005022506
31 print('datefrombufr=',datefrombufr)
32 print('filetime=',filetime)
33 #########################################################################################################
34 #########################################################################################################
35 #########################################################################################################
36 
37 hdstr1 ='SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON HOLS'
38 hdstr2 ='SAZA SOZA BEARAZ SOLAZI'
39 obstr = 'TMBR'
40 
41 ####################################################
42 # read amsua radiance file.
43 bufr = ncepbufr.open(radbufr_filename)
44 ####################################################
45 ####################################################
46 ####################################################
47 nhd1=len(hdstr1.split())
48 nhd2=len(hdstr2.split())
49 nhd=nhd1+nhd2
50 ###################################################
51 nob=0
52 nob1=0
53 
54 nrec_count=0
55 
56 # VARIABLES FROM INPUT BUFR FILE
57 DUMsaid = []
58 DUMfovn = []
59 DUMyear = []
60 DUMmnth = []
61 DUMdays = []
62 DUMhour = []
63 DUMminu = []
64 DUMseco = []
65 DUMhols = []
66 DUMclat = []
67 DUMclon = []
68 DUMsaza = []
69 DUMbearaz = []
70 DUMsoza = []
71 DUMsolazi = []
72 DUMtmbr = []
73 
74 while bufr.advance() == 0:
75  nobs_message = 0
76  print('dumsaid=', np.unique(DUMsaid))
77  while bufr.load_subset() == 0:
78  nrec_count=nrec_count+1
79  hdr1 = bufr.read_subset(hdstr1).squeeze()
80  hdr2 = bufr.read_subset(hdstr2).squeeze()
81 
82  said=int(hdr1[0])
83  fovn=int(hdr1[1])
84  year=int(hdr1[2])
85  mnth=int(hdr1[3])
86  days=int(hdr1[4])
87  hour=int(hdr1[5])
88  minu=int(hdr1[6])
89  seco=int(hdr1[7])
90  clat=float(hdr1[8])
91  clon=float(hdr1[9])
92  if clon<0:
93  clon=clon+360.
94  #print('clon=', clon)
95  hols=float(hdr1[10]) ##height of landsurface
96 
97  saza=float(hdr2[0])
98  soza=float(hdr2[1])
99  bearaz=float(hdr2[2])
100  solazi=float(hdr2[3])
101 
102  tmbr = bufr.read_subset(obstr,rep=True).squeeze()
103  nchanl = len(tmbr)
104 
105  DUMtmbr.append(tmbr)
106 
107  DUMsaid.append(said)
108  DUMfovn.append(fovn)
109  DUMyear.append(year)
110  DUMmnth.append(mnth)
111  DUMdays.append(days)
112  DUMhour.append(hour)
113  DUMminu.append(minu)
114  DUMseco.append(seco)
115  DUMhols.append(hols)
116  DUMclat.append(clat)
117  DUMclon.append(clon)
118  DUMsaza.append(saza)
119  DUMbearaz.append(bearaz)
120  DUMsoza.append(soza)
121  DUMsolazi.append(solazi)
122 
123 bufr.close()
124 print('nrec=',nrec_count)
125 #########################################################################################################
126 #########################################################################################################
127 #########################################################################################################
128 # NOTE: /home/Jeffrey.S.Whitaker/.local/lib/python2.7/site-packages/
129 # does not have satellite_names: need to use own version
130 from ncepbufr import satellite_names as code2sat_name
131 import string
132 
133 punct=string.punctuation
134 space=' '
135 hyphen='-'
136 digits=string.digits
137 
138 
139 
140 uniqsats=np.unique(DUMsaid)
141 numsat=len(uniqsats)
142 
143 for n in range(numsat):
144  holsarr = []
145  fovnarr = []
146  yeararr = []
147  mntharr = []
148  daysarr = []
149  hourarr = []
150  minuarr = []
151  secoarr = []
152  clatarr = []
153  clonarr = []
154  latarr = []
155  lonarr = []
156  clatharr = []
157  clonharr = []
158  sazaarr = []
159  bearazarr = []
160  sozaarr = []
161  solaziarr = []
162  tmbrarr = []
163 
164  satcode=uniqsats[n]
165  indcs=np.where(DUMsaid==satcode)[0]
166  print('satcode,indcs=',satcode,indcs)
167  satlite = code2sat_name[satcode]
168  print('satlite_name=',satlite)
169  satlite_lower = satlite.lower
170  sat_lower = "".join(c for c in satlite_lower())
171  satlite_nopunct = "".join(c for c in satlite_lower() if c not in punct)
172  satlite_spacepunct = "".join(c for c in satlite if c not in punct)
173  satlite_nopunctnospace = "".join(c for c in satlite_nopunct if c not in space)
174  print('lower,nopunc,spacepunc,nopuncnospace=',sat_lower,satlite_nopunct,satlite_spacepunct,satlite_nopunctnospace)
175 
176  if len(satlite_lower())>=6 and satlite_lower()[0:6]=='metop-':
177  diag_satname=satlite_lower()[0:6]
178  if satlite_lower()[6]=='2':
179  diag_satname=diag_satname + 'a'
180  if satlite_lower()[6]=='1':
181  diag_satname=diag_satname + 'b'
182  if satlite_lower()[6]=='3':
183  diag_satname=diag_satname + 'c'
184 
185  if len(satlite_nopunctnospace)>=4 and satlite_nopunctnospace[0:4]=='noaa' :
186  diag_satname='n'
187  diag_satname=diag_satname + satlite_nopunctnospace[4:]
188 
189  if len(satlite_nopunctnospace)==4 and satlite_nopunctnospace[0:4]=='aqua' :
190  diag_satname='aqua'
191 
192  print('diag_satname=', diag_satname)
193 
194  numindcs=len(indcs)
195  print('nrec_sat=',numindcs)
196 
197 
198  fovnarr=[DUMfovn[i] for i in indcs]
199  yeararr=[DUMyear[i] for i in indcs]
200  mntharr=[DUMmnth[i] for i in indcs]
201  daysarr=[DUMdays[i] for i in indcs]
202  hourarr=[DUMhour[i] for i in indcs]
203  minuarr=[DUMminu[i] for i in indcs]
204  secoarr=[DUMseco[i] for i in indcs]
205  clatarr=[DUMclat[i] for i in indcs]
206  clonarr=[DUMclon[i] for i in indcs]
207  sazaarr=[DUMsaza[i] for i in indcs]
208  bearazarr=[DUMbearaz[i] for i in indcs]
209  sozaarr=[DUMsoza[i] for i in indcs]
210  solaziarr=[DUMsolazi[i] for i in indcs]
211  holsarr=[DUMhols[i] for i in indcs]
212  tmbrarr=[DUMtmbr[i] for i in indcs]
213 
214  diagtime=filetime[0:8]+'_'+filetime[8:10]
215  netcdf_filename = netcdf_filePATH+'diag_amsua_'+diag_satname+'_anl.'+diagtime+'z.nc4'
216  print('netcdffile=',netcdf_filename)
217  nc =netCDF4.Dataset(netcdf_filename,'w',format='NETCDF4')
218  nchans_dim = nc.createDimension('nchans',None)
219  nobs_bufr_dim = nc.createDimension('nobs_bufr',None)
220  nobs_diag_dim = nc.createDimension('nobs',None)
221  Observation_Class_maxstrlen_dim = nc.createDimension('Observation_Class_maxstrlen',7)
222  BC_angord_arr_dim = nc.createDimension('BC_angord_arr_dim',5)
223 
224  chaninfoidx_nc = nc.createVariable('chaninfoidx','i','nchans') ;
225 
226  freq_nc = nc.createVariable('frequency',np.double,'nchans') ;
227 
228  polariz_nc = nc.createVariable('polarization','i','nchans') ;
229 
230  wavenum_nc = nc.createVariable('wavenumber',np.double,'nchans') ;
231  errvar_nc = nc.createVariable('error_variance',np.double,'nchans') ;
232  lapseavg_nc = nc.createVariable('mean_lapse_rate',np.double,'nchans') ;
233 
234  useflag_nc = nc.createVariable('use_flag','i','nchans') ;
235  sensorchan_nc = nc.createVariable('sensor_chan','i','nchans') ;
236  satinfochan_nc = nc.createVariable('satinfo_chan','i','nchans') ;
237  chanindx_nc = nc.createVariable('Channel_Index','i','nobs') ;
238  obclass_nc = nc.createVariable('Observation_Class','S1', ('nobs','Observation_Class_maxstrlen')) ; ###### will likelyneed netCDF4.stringtochar
239 
240 # VARIABLES FROM THE INPUT BUFR FILE (dim = nobs_bufr)
241  fovn_bufr_nc = nc.createVariable('FOVN', np.double,'nobs_bufr') ;
242  year_bufr_nc = nc.createVariable('YEAR', 'i','nobs_bufr') ;
243  mnth_bufr_nc = nc.createVariable('MNTH', 'i','nobs_bufr') ;
244  days_bufr_nc = nc.createVariable('DAYS', 'i','nobs_bufr') ;
245  hour_bufr_nc = nc.createVariable('HOUR', 'i','nobs_bufr') ;
246  minu_bufr_nc = nc.createVariable('MINU', 'i','nobs_bufr') ;
247  seco_bufr_nc = nc.createVariable('SECO', 'i','nobs_bufr') ;
248  clat_bufr_nc = nc.createVariable('CLAT', np.double,'nobs_bufr') ;
249  clon_bufr_nc = nc.createVariable('CLON', np.double,'nobs_bufr') ;
250  #clath_bufr_nc = nc.createVariable('CLATH', np.double,'nobs_bufr') ;
251  #clonh_bufr_nc = nc.createVariable('CLONH', np.double,'nobs_bufr') ;
252  hols_bufr_nc = nc.createVariable('HOLS', np.double,'nobs_bufr') ;
253 
254  satzen_bufr_nc = nc.createVariable('SAZA', np.double,'nobs_bufr') ;
255  solzen_bufr_nc = nc.createVariable('SOZA', np.double,'nobs_bufr') ;
256  bearaz_bufr_nc = nc.createVariable('BEARAZ', np.double,'nobs_bufr') ;
257  solaz_bufr_nc = nc.createVariable('SOLAZI', np.double,'nobs_bufr') ;
258 
259  tmbr_bufr_nc= nc.createVariable('TMBR', np.double,('nobs_bufr','nchans')) ;
260 
261 # VARIABLES FROM THE OBS DIAG FILE (dim = nobs)
262 
263 
264  lat_diag_nc = nc.createVariable('Latitude', np.double,'nobs') ;
265  lon_diag_nc = nc.createVariable('Longitude', np.double,'nobs') ;
266  elev_diag_nc = nc.createVariable('Elevation', np.double,'nobs') ;
267  time_diag_nc = nc.createVariable('Obs_Time', np.double,'nobs') ;
268  scanpos_diag_nc = nc.createVariable('Scan_Position', np.double,'nobs') ;
269  satzen_diag_nc = nc.createVariable('Sat_Zenith_Angle', np.double,'nobs') ;
270  satazi_diag_nc = nc.createVariable('Sat_Azimuth_Angle', np.double,'nobs') ;
271  solzen_diag_nc = nc.createVariable('Sol_Zenith_Angle', np.double,'nobs') ;
272  solazi_diag_nc = nc.createVariable('Sol_Azimuth_Angle', np.double,'nobs') ;
273  glintang_diag_nc = nc.createVariable('Sun_Glint_Angle', np.double,'nobs') ;
274  waterfrac_diag_nc= nc.createVariable('Water_Fraction', np.double,'nobs') ;
275  landfrac_diag_nc = nc.createVariable('Land_Fraction', np.double,'nobs') ;
276  icefrac_diag_nc = nc.createVariable('Ice_Fraction', np.double,'nobs') ;
277  snowfrac_diag_nc = nc.createVariable('Snow_Fraction', np.double,'nobs') ;
278  waterT_diag_nc = nc.createVariable('Water_Temperature', np.double,'nobs') ;
279  landT_diag_nc = nc.createVariable('Land_Temperature', np.double,'nobs') ;
280  iceT_diag_nc = nc.createVariable('Ice_Temperature', np.double,'nobs') ;
281  snowT_diag_nc = nc.createVariable('Snow_Temperature', np.double,'nobs') ;
282  soilT_diag_nc = nc.createVariable('Soil_Temperature', np.double,'nobs') ;
283  soilM_diag_nc = nc.createVariable('Soil_Moisture', np.double,'nobs') ;
284  #### INT ####
285  landTYPE_diag_nc = nc.createVariable('Land_Type_Index','i','nobs') ;
286  #### INT ####
287  tsavg5_diag_nc= nc.createVariable('tsavg5', np.double,'nobs') ;
288  sstcu_diag_nc= nc.createVariable('sstcu', np.double,'nobs') ;
289  sstph_diag_nc= nc.createVariable('sstph', np.double,'nobs') ;
290  sstnv_diag_nc= nc.createVariable('sstnv', np.double,'nobs') ;
291  dta_diag_nc= nc.createVariable('dta', np.double,'nobs') ;
292  dqa_diag_nc= nc.createVariable('dqa', np.double,'nobs') ;
293  dtp_avh_diag_nc= nc.createVariable('dtp_avh', np.double,'nobs') ;
294  vegfrac_diag_nc= nc.createVariable('Vegetation_Fraction', np.double,'nobs') ;
295  snodep_diag_nc= nc.createVariable('Snow_Depth', np.double,'nobs') ;
296  tpwc_diag_nc= nc.createVariable('tpwc_amsua', np.double,'nobs') ;
297  clwges_diag_nc= nc.createVariable('clw_guess_retrieval', np.double,'nobs') ;
298  sfcwind_diag_nc= nc.createVariable('Sfc_Wind_Speed', np.double,'nobs') ;
299  cloudfrac_diag_nc= nc.createVariable('Cloud_Frac', np.double,'nobs') ;
300  CTP_diag_nc= nc.createVariable('CTP', np.double,'nobs') ;
301  CLW_diag_nc= nc.createVariable('CLW', np.double,'nobs') ;
302  TPWC_diag_nc= nc.createVariable('TPWC', np.double,'nobs') ;
303  clwobs_diag_nc= nc.createVariable('clw_obs', np.double,'nobs') ;
304  clwges_diag_nc= nc.createVariable('clw_guess', np.double,'nobs') ;
305  foundT_diag_nc= nc.createVariable('Foundation_Temperature', np.double,'nobs') ;
306  sst_warmdt_diag_nc= nc.createVariable('SST_Warm_layer_dt', np.double,'nobs') ;
307  sst_cooldt_diag_nc= nc.createVariable('SST_Cool_layer_tdrop', np.double,'nobs') ;
308  sst_dtzdtf_diag_nc= nc.createVariable('SST_dTz_dTfound', np.double,'nobs') ;
309  obs_diag_nc= nc.createVariable('Observation', np.double,'nobs') ;
310  oMf_adj_diag_nc= nc.createVariable('Obs_Minus_Forecast_adjusted', np.double,'nobs') ;
311  oMf_unadj_diag_nc= nc.createVariable('Obs_Minus_Forecast_unadjusted', np.double,'nobs') ;
312  invOBSerr_diag_nc= nc.createVariable('Inverse_Observation_Error', np.double,'nobs') ;
313  qcflag_diag_nc= nc.createVariable('QC_Flag', np.double,'nobs') ;
314  emiss_diag_nc= nc.createVariable('Emissivity', np.double,'nobs') ;
315  wtd_lapse_diag_nc= nc.createVariable('Weighted_Lapse_Rate', np.double,'nobs') ;
316  dtb_dts_diag_nc= nc.createVariable('dTb_dTs', np.double,'nobs') ;
317  bc_const_diag_nc= nc.createVariable('BC_Constant', np.double,'nobs') ;
318  bc_scanangle_diag_nc= nc.createVariable('BC_Scan_Angle', np.double,'nobs') ;
319  bc_clw_diag_nc= nc.createVariable('BC_Cloud_Liquid_Water', np.double,'nobs') ;
320  bc_lapsesqd_diag_nc= nc.createVariable('BC_Lapse_Rate_Squared', np.double,'nobs') ;
321  bc_lapse_diag_nc= nc.createVariable('BC_Lapse_Rate', np.double,'nobs') ;
322  bc_cosinelatxnode_diag_nc= nc.createVariable('BC_Cosine_Latitude_times_Node', np.double,'nobs') ;
323  bc_sine_diag_nc= nc.createVariable('BC_Sine_Latitude', np.double,'nobs') ;
324  bc_emiss_diag_nc= nc.createVariable('BC_Emissivity', np.double,'nobs') ;
325  bc_fixedscanpos_diag_nc= nc.createVariable('BC_Fixed_Scan_Position', np.double,'nobs') ;
326  bc_angord_diag_nc= nc.createVariable('BC_angord', np.double,('nobs','BC_angord_arr_dim')) ;
327 
328  #// global attributes:
329  ##nc.description='blahblah'
330  nc.Satellite_Sensor = 'amsua_'+diag_satname ;
331  nc.Satellite = diag_satname ;
332  nc.Observation_type = "amsua" ;
333 # should not be hard coded
334  nc.Outer_Loop_Iteration = -1;
335  nc.Number_of_channels = -1 ;
336  nc.Number_of_Predictors = -1 ;
337  nc.date_time = -1 ;
338  nc.ireal_radiag = -1 ;
339  nc.ipchan_radiag = -1 ;
340  nc.iextra = -1 ;
341  nc.jextra = -1 ;
342  nc.idiag = -1 ;
343  nc.angord = -1 ;
344  nc.iversion_radiag = -1 ;
345  nc.New_pc4pred = -1 ;
346  nc.ioff0 = -1 ;
347 
348 
349 
350  #########################################################################################################
351  fovnarr = np.array(fovnarr)
352  clatarr = np.array(clatarr)
353  clonarr = np.array(clonarr)
354  latarr = -9999*np.ones(clatarr.shape)
355  lonarr = -9999*np.ones(clonarr.shape)
356  clatharr = -9999*np.ones(clatarr.shape)
357  clonharr = -9999*np.ones(clonarr.shape)
358  holsarr = np.array(holsarr)
359 
360  sazaarr = np.array(sazaarr)
361  bearazarr = np.array(bearazarr)
362  sozaarr = np.array(sozaarr)
363  solaziarr = np.array(solaziarr)
364  tmbrarr = np.array(tmbrarr)
365  #print('clatarray=',clatarr)
366 
367  nc['FOVN'][:] = fovnarr
368  nc['YEAR'][:] = yeararr
369  nc['MNTH'][:] = mntharr
370  nc['DAYS'][:] = daysarr
371  nc['HOUR'][:] = hourarr
372  nc['MINU'][:] = minuarr
373  nc['SECO'][:] = secoarr
374  nc['CLAT'][:] = clatarr
375  nc['CLON'][:] = clonarr
376  nc['HOLS'][:] = holsarr
377  #nc['CLATH'][:] = clatharr
378  #nc['CLONH'][:] = clonharr
379  nc['SAZA'][:] = sazaarr
380  nc['SOZA'][:] = sozaarr
381  nc['BEARAZ'][:] = bearazarr
382  nc['SOLAZI'][:] = solaziarr
383 
384  nc['TMBR'][:] = tmbrarr
385  nc.sync() # dump data to disk.
386 
387 
388  nc.close()
389  del nc
390  del indcs
391  del fovnarr
392  del latarr
393  del lonarr
394  del clatarr
395  del clonarr
396  del clatharr
397  del clonharr
398  del sazaarr
399  del bearazarr
400  del sozaarr
401  del solaziarr
402  del holsarr
403  del tmbrarr
404 
405 
406 
407 
408 
409 
410 
411 
JEDI_amsu_bufr2nc.satlite_lower
satlite_lower
Definition: JEDI_amsu_bufr2nc.v8.py:169