GRIB2 ファイルを自力でパースしてみた [python]

京都大学生存圏研究所 生存圏データベース > グローバル大気観測データ > 気象庁データ > 全国合成レーダーGPV のデータを Python 上で扱ってみます。

この全国合成レーダーGPVは「気象庁が保有する全国20台の気象レーダーで観測したエコー強度(レーダーで観測される換算降水強度)とエコー頂高度(レーダーで観測される降水エコーの高さ)」で10分毎に提供されています。

ファイルの形式は GRIB2 です。
GRIB2（国際気象通報式 FM92 GRIB 二進形式格子点資料気象通報式（第2版）） は気象の格子データの格納に使われるフォーマットです。

Python で GRIB2 を読み出すライブラリは ecCodes など様々あるようですが，うまく扱えなかったので今回は自力でパースしようと思います。
そして，手順を備忘録がてら残しておきます。

データ利用上の注意点

京都大学生存圏研究所 生存圏データベース より

• 本データベースのデータの使用にあたっては、それぞれのデータベースのページにある注意事項をお守りください。
• 本データベースの作成には、日本学術振興会科学研究費補助金「研究成果公開促進費」が利用されています。
• 出版物（科学論文を含む）やプロダクトに利用する場合には、 次の文章を引用していただきますよう、よろしくお願いいたします。
  「これらのデータは京都大学生存圏研究所が運営する生存圏データベースによって収集・配布されたものです（http://database.rish.kyoto-u.ac.jp）。」
  （英語の場合）
  "These data were collected and distributed by Research Institute for Sustainable Humanosphere, Kyoto University (http://database.rish.kyoto-u.ac.jp/index-e.html)."
  本データベースから取得したデータを利用した科学論文が出版された際には、著者名、タイトル、論文誌、巻号、ページなどを ***** までお送りください。
  データベースの維持・発展のため、ご協力をお願いいたします。

気象庁データ より

• (財)気象業務支援センターを通して公開されている気象庁作成の数値予報データ・観測データのダウンロードサイトです．
• ここでは教育研究機関向けにデータを提供しています．企業活動等のためにデータを頻繁に必要とされる方は，気象業務支援センターからデータを直接購入し，データ提供スキーム全体の維持発展にご協力ください．
• 本データの収集は KAGI21 の活動の一環として行われました．

教育研究機関向けだそうです。
個人的に遊びたいときも購入しなきゃいけないんでしょうか…？
まぁ，それは置いておいて次へ進みます。

自力パースの前準備

フォーマットの詳細は 配信資料に関する技術情報（気象編）第 162 号 ～1km メッシュ全国合成レーダーGPV の提供等について～ にあります。
また，データ中のランレングス圧縮については ランレングス符号化法の解説 に解説があります。

これらの資料を参考に実際に 1km メッシュのデータを読み出すスクリプトを用意します。

スクリプト

from datetime import datetime
import numpy as np

class JMAGPV(object):
    def __init__(self, file):
        self.f = file
        self.read()

    def read(self):
        pass
        
    def r(self, start, end=-1):
        assert self.f != None
        if end == -1:
            end = start
        return self.f.read(end-start+1)
        
    def r_int(self, start, end=-1, byte=None):
        return int.from_bytes(self.r(start, end) if byte == None else byte, 'big')

    def r_str(self, start, end=-1):
        return self.r(start, end).decode()
    
    def r_dt(self, start):
        return datetime(self.r_int(start, start+1),
                        self.r_int(start+2),
                        self.r_int(start+3),
                        self.r_int(start+4),
                        self.r_int(start+5),
                        self.r_int(start+6))
    @staticmethod
    def rle_repeat(rle, lngu, maxv):
        return sum([(lngu**(m))*(n - (maxv+1)) for m, n in enumerate(rle)])+1 if rle != [] else 1
    
    @staticmethod
    def decode_rle(b, nbit, maxv, column=2560, row=3360, start_row=0, end_row=3360, start_column=0, end_column=2560):
        prev = None
        rle = []
        result = []
        lngu = 2**nbit - 1 - maxv
        for i in b:
            if i <= maxv:
                if prev != None:
                    result.extend([prev] * JMAGPV.rle_repeat(rle, lngu, maxv))
                    rle = []
                prev = i
            else:
                rle += [i]
        result.extend([prev] * JMAGPV.rle_repeat(rle, lngu, maxv))
        return np.array(result, dtype=np.uint8).reshape(row, column)

class Ggis1km(JMAGPV):
    class Section0:
        name = 0
        def __init__(self, grib, _, document_field, grib_version, bytes_length):
            self.grib = grib
            self.document_field = document_field
            self.grib_version = grib_version
            self.bytes_length = bytes_length

    class Section1:
        name = 1
        def __init__(self, length, center, secondary_center, grib_master_version, grib_local_version,
                     reference_time, document_datetime, create_status, document_type):
            self.length = length
            self.center = center
            self.secondary_center = secondary_center
            self.grib_master_version = grib_master_version
            self.grib_local_version = grib_local_version
            self.reference_time = reference_time
            self.document_datetime = document_datetime
            self.create_status = create_status
            self.document_type = document_type

    class Section3:
        name = 3
        def __init__(self, length, lattice_system, num_of_points, _a, _b, lattice_system_definition_template_num,
                     earth_shape, _c, _d,
                     long_axis_scale_factor, long_axis_scale_length, short_axis_scale_factor, short_axis_scale_length,
                     lat_points, lng_points, base_angle, _e, first_lat, first_lng, resolution_and_component_flag,
                     last_lat, last_lng, incremental_i, incremental_j, scanning_mode):
            self.length = length
            self.lattice_system = lattice_system
            self.num_of_points = num_of_points
            self.lattice_system_definition_template_num = lattice_system_definition_template_num
            self.earth_shape = earth_shape
            self.long_axis_scale_factor = long_axis_scale_factor
            self.long_axis_scale_length = long_axis_scale_length
            self.short_axis_scale_factor = short_axis_scale_factor
            self.short_axis_scale_length = short_axis_scale_length
            self.lat_points = lat_points
            self.lng_points = lng_points
            self.base_angle = base_angle
            self.first_lat = first_lat
            self.first_lng = first_lng
            self.resolution_and_component_flag = resolution_and_component_flag
            self.last_lat = last_lat
            self.last_lng = last_lng
            self.incremental_i = incremental_i
            self.incremental_j = incremental_j
            self.scanning_mode = scanning_mode

    class Section4:
        name = 4
        def __init__(self, length, _a, product_template_num, param_category, param_num, processing_type,
                     value_type, _b, _c, _d, _e, _f, first_fixed_surface_type, _g, _h, _i, _j, _k,
                     all_timedelta_finish_datetime, _l, _m, statistical_processing_type,
                     statistical_processing_unit, statistical_processing_timedelta, _n, _o, _p,
                     radar_operation_info, rf_conversion_factor_operation_info, _q):
            self.length = length
            self.product_template_num = product_template_num
            self.param_category = param_category
            self.param_num = param_num
            self.processing_type = processing_type
            self.value_type = value_type
            self.first_fixed_surface_type = first_fixed_surface_type
            self.all_timedelta_finish_datetime = all_timedelta_finish_datetime
            self.statistical_processing_type = statistical_processing_type
            self.statistical_processing_unit = statistical_processing_unit
            self.statistical_processing_timedelta = statistical_processing_timedelta
            self.radar_operation_info = radar_operation_info
            self.rf_conversion_factor_operation_info = rf_conversion_factor_operation_info

    class Section5:
        name = 5
        def __init__(self, length, num_of_all_document_points, document_template_num, num_of_bits_by_one_data,
                     max_level_value, _a, data_representative_value, level_data_representative_value):
            self.length = length
            self.num_of_all_document_points = num_of_all_document_points
            self.document_template_num = document_template_num
            self.num_of_bits_by_one_data = num_of_bits_by_one_data
            self.max_level_value = max_level_value
            self.data_representative_value = data_representative_value
            self.level_data_representative_value = level_data_representative_value

    class Section6:
        name = 6
        def __init__(self, length, bitmap_indicator):
            self.length = length
            self.bitmap_indicator = bitmap_indicator

    class Section7:
        name = 7
        def __init__(self, length, run_length_octet_strings):
            self.length = length
            self.run_length_octet_strings = run_length_octet_strings

    def read(self):
        b  = self.r
        i  = self.r_int
        s  = self.r_str
        dt = self.r_dt

        # 第 0 節
        self.section0 = Ggis1km.Section0(s(1,4), i(5,6), i(7), i(8), i(9, 16))
        # 第 1 節
        sec_length = i(1,4)
        sec_name = i(5)
        if sec_name == 1:
            self.section1 = Ggis1km.Section1(sec_length, i(6,7), i(8,9), i(10), i(11), i(12), dt(13), i(20), i(21))
        # 第 3 節
        sec_length = i(1,4)
        sec_name = i(5)
        if sec_name == 3:
            self.section3 = Ggis1km.Section3(sec_length, i(6), i(7, 10), i(11), i(12), i(13,14), i(15), i(16),
                                             i(17, 20), i(21),i(22, 25), i(26), *[i(27, 30) for _ in range(7)],
                                             i(55), *[i(56, 59) for _ in range(4)], i(72))
        # 第 4 節
        sec_length = i(1,4)
        sec_name = i(5)
        if sec_name == 4:
            self.section4 = Ggis1km.Section4(sec_length, i(6,7), i(8,9), *[i(10) for _ in range(10, 15)],
                                             i(15,16), i(17), i(18), i(19,22), i(23), i(24), i(25,28), i(29), i(30),
                                             i(31,34), dt(35), i(42), i(43,46), i(47), i(48), i(49), i(50,53), i(54),
                                             i(55,58), *[i(59,66) for _ in range(3)])
        # 第 5 節
        sec_length = i(1,4)
        sec_name = i(5)
        if sec_name == 5:
            self.section5 = Ggis1km.Section5(sec_length, i(6,9), i(10,11), i(12), i(13,14), i(15,16), i(17), 
                                             [i(16+2*nn,17+2*nn) for nn in range(1, 252)])
        # 第 6 節
        sec_length = i(1,4)
        sec_name = i(5)
        if sec_name == 6:
            self.section6 = Ggis1km.Section6(sec_length, i(6))
        # 第 7 節
        sec_length = i(1,4)
        sec_name = i(5)
        if sec_name == 7:
            self.section7 = Ggis1km.Section7(sec_length, b(6,sec_length))
        self.data = self.decode_rle(self.section7.run_length_octet_strings, 8, self.section5.max_level_value)
        flag = b(1,4)
        assert flag.decode() == '7777'

from pathlib import Path
with Path(Path()/'data'/'Z__C_RJTD_20070101000000_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.bin').open('rb') as f:
    g = Ggis1km(f)

で読み出せるようになりました（多分）。
時間がかかりすぎて検証できてないですが。

スクリプトのミスを修正したら無事動きました！

気象庁 テンプレート 7.200 ランレングス圧縮の解凍

def rle_repeat(rle, lngu, maxv):
    return sum([(lngu**(m))*(n - (maxv+1)) for m, n in enumerate(rle)])+1 if rle != [] else 1

def decode_rle(b, nbit, maxv, column=2560, row=3360, start_row=0, end_row=3360, start_column=0, end_column=2560):
    prev = None
    rle = []
    result = []
    lngu = 2**nbit - 1 - maxv
    for i in b:
        if i <= maxv:
            if prev != None:
                result.extend([prev] * rle_repeat(rle, lngu, maxv))
                rle = []
            prev = i
        else:
            rle += [i]
    result.extend([prev] * rle_repeat(rle, lngu, maxv))
    return np.array(result, dtype=np.uint8).reshape(row, column)