Postprocessing of libecbufr to make stats of TEMP BUFRS

posttemp.rb

#!/usr/bin/ruby

class BufrMsg

  def initialize
    @head = {}
    @edition = nil
  end

  def inject key, val
    case val
    when /^"/ then
      case val
      when /([A-Z]{4}\d{0,2}\\040[A-Z]{4})\\040(\d{6})(\\040[A-Z]{3})?/
      #    /1............................1.....2.....23.............3./
	taic, ygg, bbb = $1, $2, $3
	@head['taic'] = taic.sub(/\\040/, '')
	@head['ygg'] = ygg
	@head['bbb'] = bbb ? bbb.sub(/\\040/, '') : '///'
      when /Z__C_RJTD_(\d{14})_/
	t = $1
	@head['taic'] = "file__RJTD"
	@head['ygg'] = t[6,6]
	@head['bbb'] = '///'
      else
        puts "# #{key}=#{val}" if $VERBOSE
      end
    when /^\d+$/
      @head[key] = val.to_i
    end
  end

  def compile
    @edition = @head['BUFR_EDITION']
    t = @head.values_at(*%w(YEAR MONTH DAY HOUR MINUTE SECOND))
    t[0] += 2000 if @edition <= 3
    @head['time'] = Time.gm(*t)
    @head['stime'] = @head['time'].strftime('%Y%m%dT%H%M%S')
    c, s = @head['ORIG_CENTER'], @head['ORIG_SUB_CENTER']
    s = 65535 if s.nil?
    @head['center'] = [
      c == 65535 ? '///' : format('%03u', c),
      s == 65535 ? '///' : format('%03u', s)
    ].join
    @head['usn'] = format('%02u', @head['UPDATE_SEQUENCE'])
    c = @head.values_at(*%w(DATA_CATEGORY
      INTERN_SUB_CATEGORY LOCAL_SUB_CATEGORY))
    @head['cat'] = if @edition >= 4
      then format('%02X%02X%02X', *c)
      else format('%02X//%02X', c[0], c[2])
      end
    c = @head.values_at(*%w(MASTER_TABLE_VERSION
      LOCAL_TABLE_VERSION))
    @head['tbv'] = format('%02X%02X', *c)
    @head.freeze
  rescue Interrupt => e
    STDERR.puts @head.inspect
    raise e
  end

  def report
    format('%-10s,%6s,%3s,Ed%1u,Tv%04s,Cat%06s,Ctr%06s,%15s,Usn%02u',
      @head['taic'], @head['ygg'], @head['bbb'],
      @edition, @head['tbv'],
      @head['cat'], @head['center'], @head['stime'], @head['usn']
      )
  end

  def sstime
    @head['time'].strftime('%dT%H')
  end

  def taic
    @head['taic']
  end

  def edition
    @edition
  end

  def center
    @head['center']
  end

end

class BufrSubset

  def initialize msg
    @msg = msg
    @d1 = @d2 = nil
    @dcd = {}
    @hack3 = @hackr = nil
    @prevr = nil
    @diag = {}
  end

  def taic
    @msg.taic
  end

  def inject_main desc, r, data
    #printf ": %06s %-5s %s\n", desc, r, data if $DEBUG
    if /^3/ =~ desc and (/^031/ =~ @d1 and /^1..000/ =~ @d2 or /^1/ =~ @d1 or /^102/ =~ @d2)
      @hack3 = desc
      @hackr = "0"
      #STDERR.puts "hack trigger #{desc}" if $DEBUG
    end
    if desc == @hack3
      @hackr = @hackr.next
      #STDERR.puts "hack increment #{@hackr}" if $DEBUG
    end
    if @hackr and r == '0'
      r = @hackr
    end
    dcd = @dcd
    if r
      n = 0
      r.split(/\./).map{|s|
        n <<= 14
	n |= s.to_i
      }
      if r != @prevr
        dcd[n | (1 << 27)] = dcd[n] if dcd[n]
	dcd[n] = {}
        #STDERR.puts "alloc #{n}" if $DEBUG 
      end
      dcd = dcd[n]
      @prevr = r
    end
    dcd[desc] = data
  end

DPAT = /^(?#1=)([0-3][0-9]{5})(?: \{R=(?#2=)([.0-9]+)\}(?:\{[^}]*\})?)? (?:\((?#3=)(0x\w+):\d+bits\))?(?#4=)("[^"]*"|[^"]\S*)?/

  def inject line
    case line
    when /^#([0-3][0-9]{5}) $/ then
      desc = $1
      #@dcd[desc] = false
    when DPAT
      desc, r, abit, data = $1, $2, $3, $4
      if data
        data = 'MSNG' if abit and abit.to_i != 0
	data = case data
	  when /^"([^"]*)"/ then $1
	  when 'MSNG' then nil
	  when /\./ then data.to_f
	  when /^[01]+$/ 
	    if /^008(001|042)/ =~ desc
	    then data.to_i(2)
	    else data.to_i
	    end
	  else data.to_i
	  end
      end
      inject_main(desc, r, data)
      @d2 = @d1
      @d1 = desc
    else
      STDERR.puts "unrecognised <#{line}>"
      return nil
    end
  end

  def each_lev
    @dcd.keys.select{|i| Integer === i}.sort.each {|i|
      yield @dcd[i]
    }
  end

  def each_p
    pprev = 0
    each_lev { |v|
      next unless p = v['007004']
      p = p.to_i
      if p < 50_00 and pprev > 1000_00 and /^034000/ === @msg.center
        STDERR.puts "# RJTD-P special #{p}"
        p += 1000_00
      end
      yield v, p
      pprev = p
    }
  end

  def vslice desc
    r = []
    each_lev {|v|
      r.push v[desc] if v.include?(desc)
    }
    r
  end

  STDP_A = [1000_00, 925_00, 850_00, 700_00, 500_00,
    400_00, 300_00, 250_00, 200_00, 150_00]
    # 100_00 omitted intentionally
  STDP_C = [70_00, 50_00, 30_00, 20_00, 10_00,
    7_00, 5_00, 3_00, 2_00, 1_00]
  STDP = {}
  (STDP_A + STDP_C).each{|p| STDP[p] = true}

  def diag_pval
    pres = []
    each_p {|v, p| pres.push p }
    @diag['np'] = pres.size
    @diag['pmax'] = begin pres.compact.max / 100 rescue nil end
    @diag['pmin'] = begin pres.compact.min / 100 rescue nil end
  end

  def diag_zval
    hgt = vslice('007003')
    @diag['nz'] = hgt.size
    @diag['zmin'] = begin hgt.compact.min rescue '#####' end
    @diag['zmax'] = begin hgt.compact.max rescue '#####' end
  end

  def diag_part2
    nt = nw = no = nx = 0
    ps = []
    # p loop
    each_p {|v, p|
      t = (v['012001'] || v['012101'])
      d = v['011001']
      if d and d > 360
        STDERR.puts "# RPMM special dd=#{d}" # dd=381 for unknown reason
        d = nil
      end
      vs = (v.include?('011061') || v.include?('011062'))
      if STDP[p] and (v.include?('010009') || v.include?('010008') || v.include?('010003')) then
	ps.push p
      elsif vs then
        nx += 1
      elsif t and not d then # t sig levs
        nt += 1
      elsif d and not t then # wind levels incl. maxw
        nw += 1
      else
        no += 1
      end
    }
    if ps.size > 0 and nx <= 3 and no <= 4 and nt + nw < 7 then
      if ps.min >= 100_00 then
        @diag['P'] = 'TA'
      elsif ps.max <= 100_00 then
        @diag['P'] = 'TC'
      else
        @diag['P'] = 'TAC'
      end
    elsif ps.size <= 3 then
      if @diag['pmin'] >= 100 then
        @diag['P'] = 'TB'
      elsif @diag['pmax'] <= 100 then
        @diag['P'] = 'TD'
      else
        @diag['P'] = 'TBD'
      end
    else
      if @diag['pmin'] >= 100 then
        @diag['P'] = 'TAB'
      elsif @diag['pmax'] <= 100 then
        @diag['P'] = 'TCD'
      else
        @diag['P'] = 'TABCD'
      end
    end
    STDERR.puts "# #{@diag['P']} ps #{ps.inspect} nx #{nx} no #{no} nt #{nt} nw #{nw}" if $VERBOSE
  end

  def diag_part1
    tstat = Hash.new(0)
    nzhi = nzlo = 0
    ndhi = ndlo = 0
    nmov = 0
    # p loop
    each_p {|v, p|
      if t = (v['012001'] || v['012101']) then
        tsub = (t * 100 + 0.5).floor
	tstat[tsub % 20] += 1
      end
      if p <= 500_00 then
        z = v['010009'] || v['010008'] || v['010003']
	if /^028/ === @msg.center and p < 200_00 and z and z < 1000
	  STDERR.puts "# DEMS-nzhi special"
	  z = nil
	end
	if z.nil? then
	  :skip
        elsif (z % 10) == 0 then
	  nzlo += 1
	else
	  nzhi += 1
	end
      end
      if d = v['011001'] and d < 360
        if (d % 5) == 0 then
	  ndlo += 1
	else
	  #STDERR.puts v.inspect
	  ndhi += 1
	end
      end
      # DEMS fills constant zero in those fields in reformatted bufr.
      if v['004086'].to_f != 0.0 and v['005015'].to_f != 0.0 and v['006015'].to_f != 0.0 then
        nmov += 1
      end
    }
    tpat = case tstat.size
      when 0 then nil
      when 1 then 1
      when 2 then 1
      when 3 then
        STDERR.puts "# DEMS-tpat special"
        1
      else 20
      end
    if @diag['zmin'] and @diag['zmax'] then
      if @diag['zmin'] >= 15000 then
        @diag['P'] = 'ZH'
      elsif @diag['zmax'] < 19000 then
        @diag['P'] = 'ZL'
      else
        @diag['P'] = 'ZLH'
      end
    elsif tpat == 20 or nzhi > 0 or ndhi > 1 then
      STDERR.puts "# too short Native BUFR - nt #{tstat.inspect} nz #{nzhi}/#{nzlo} nd #{ndhi}/#{ndlo}" if @diag['np'] < 150 and $VERBOSE
      if @diag['pmax'] < 100 then
        @diag['P'] = 'NH'
      elsif @diag['pmin'] >= 100 then
        @diag['P'] = 'NL'
      else
        @diag['P'] = 'NLH'
      end
    elsif tpat == 1 or nzlo > 0 or ndlo > 0 then
      diag_part2
    else
      @diag['P'] = '/'
    end
    if nmov * 2 > ndhi then
      @diag['P'] = ['+', @diag['P']].join
    end
  end

  def diag_sgnf
    sfc = std = trp = mxw = sgt = sgw = uid = 0
    each_lev {|v|
      if s = v["008001"]
        sfc += 1 unless (s & 0x40).zero?
        std += 1 unless (s & 0x20).zero?
        trp += 1 unless (s & 0x10).zero?
        mxw += 1 unless (s & 0x08).zero?
        sgt += 1 unless (s & 0x04).zero?
        sgw += 1 unless (s & 0x02).zero?
        uid += 1 if (s & 0x7E).zero?
      elsif s = v["008042"]
        sfc += 1 unless (s & 0x20000).zero?
        std += 1 unless (s & 0x10000).zero?
        trp += 1 unless (s & 0x08000).zero?
        mxw += 1 unless (s & 0x04000).zero?
        sgt += 1 unless (s & 0x03000).zero?
        sgw += 1 unless (s & 0x00800).zero?
        uid += 1 if (s & 0x3F800).zero?
      end
    }
    @diag['nsfc'] = sfc
    @diag['nstd'] = std
    @diag['ntrp'] = trp
    @diag['nmxw'] = mxw
    @diag['nsgt'] = sgt
    @diag['nsgt'] = sgt
    @diag['nnil'] = uid
  end

  def diag_stid
    begin
      i = @dcd['001001'] * 1000 + @dcd['001002']
      i = format('%05u', i)
    rescue NoMethodError
      i = nil
    end
    if i.nil? and n = @dcd['001011']
      case n
      when /^[A-Z][\x20-\x7E]*$/ then i = n.strip
      end
    end
    if i.nil?
      la, lo = @dcd['005001'].to_i, @dcd['006001'].to_i
      i = [
        format('%02u', la.abs),
	la >= 0 ? 'N' : 'S',
        format('%03u', lo.abs),
	lo >= 0 ? 'E' : 'W'
      ].join
    end
    @diag['Stid'] = i
  end

  def useless?
    return true if @diag['Stid'] == '00N000E'
    return true if @diag['P'] == '/'
    false
  end

  def compile
    p @dcd if $DEBUG
    diag_stid
    diag_pval
    diag_zval
    diag_part1
    #diag_sgnf
  end

  def stid
    @diag['Stid']
  end

  def reptype
    [@diag['P'], @msg.edition, '(', @msg.taic, ')'].join
  end

  def ptype
    [@diag['P'], @msg.edition].join
  end

  def myreport
    r = @diag.keys.sort.map{|k|
      case v = @diag[k]
      when Integer then v.to_s
      else v.inspect
      end
      [k, v].join('=')
    }
    r.unshift('')
    r.push('')
    r.join(";")
  end

  def report
    [@msg.report, myreport].join(',')
  end

  def self.parse fp
    state = :init
    line = msg = subset = nil
    for line in fp
      line.chomp!
      case state
      when :init
	case line
	when /^Error/ then :ignore
	when /^(\w+)=/
	  msg = BufrMsg.new
	  msg.inject($1, $')
	  state = :head
	when /^DATASUBSET \d+/
	  subset = BufrSubset.new(msg)
	  state = :main
	end
      when :head
	case line
	when /^(\w+)=/ then msg.inject($1, $')
	when /^DATASUBSET \d+/
	  msg.compile
	  subset = BufrSubset.new(msg)
	  state = :main
	end
      when :main
	case line
	when /^$/
	  subset.compile
	  yield subset
	  state = :init
	else
	  subset.inject(line)
	end
      end
    end
  rescue Interrupt => e
    puts line
    raise e
  end

  def self.dump fp
    db = {}
    rs = Hash.new(0)
    hs = {}
    vs = {}
    begin
    parse(fp) {|subset|
      next if subset.useless?
      STDERR.puts subset.report if $VERBOSE
      stid = subset.stid
      db[stid] = {} unless db[stid]
      db[stid][subset.reptype] = true
      rs[subset.ptype] += 1
      hs[subset.taic] = {} unless hs[subset.taic]
      hs[subset.taic][subset.ptype] = true
      vs[subset.taic] = {} unless vs[subset.taic]
      vs[subset.taic][stid] = true
    }
    rescue Interrupt
    end
    puts "\f#,parttypes"
    for ptype in rs.keys.sort
      puts [ptype, rs[ptype]].join(',')
    end
    puts "\f#,parttypes by headers"
    for taic in hs.keys.sort
      puts [taic,  ':,', hs[taic].keys.sort.join(',')].join
    end
    puts "\f#,stations by headers"
    for taic in vs.keys.sort
      puts [taic,  ':,', vs[taic].keys.sort.join(',')].join
    end
    puts "\f#,parttypes by stations"
    for stid in db.keys.sort
      puts [stid, ':,', db[stid].keys.sort.join(',')].join
    end
  end

end

ARGF.binmode
BufrSubset.dump(ARGF) if $0 == __FILE__