<jupyter_start><jupyter_text>Digging into GRIB2 binary files (JMA landslide risk)*Important:* In this repository, we provide some scripts to handle the numerical landslide risk assessment figures prepared in 10-minute intervals by the Japan Meteorological Agency (JMA). This data is  *not freely available data*, but must rather be purchased from the Japan Meteorological Business Support Center (JMBSC), which operates under the oversight of the JMA. Here are a couple relevant links:- Up to date visualized plot of landslide risk figures: https://www.jma.go.jp/jp/doshamesh/- JMBSC page for landslide risk data set: http://www.jmbsc.or.jp/jp/online/file/f-online60210.htmlThis data is only provided in a standard binary format used frequently in meteorology, called *FM92 GRIB edition 2*. The software in this repository assumes that the user has already acquired such data via the appropriate channels.__Contents:__- Opening and reading the GRIB file- Section 0: 指示節 (instructions)- Section 1: 識別節 (technical specifications)- Section 2: 地域使用節 (region specifications)- Section 3: 格子系定義節 (lattice definitions)- Section 4: プロダクト定義節 (product definitions)- Section 5: 資料表現節 (formatting)- Section 6: ビットマップ節 (bitmap info)- Section 7: 資料節 (data)  - Aside: JMA run-length compression algorithm  - Decompression of lattice data- Section 8: 終端節 (closing section)___  Opening and reading the GRIB fileTo begin, let's clarify two points: (a) the nature of the data we are looking at; (b) the format of this data.Regarding the nature of the data of interest, we are focusing exclusively on the landslide risk information computed and published by the JMA. The Japanese name is 土砂災害警戒判定メッシュ情報 (see links at the top of this page). The JMA assigns a numerical score to each point on a lattice, corresponding to the risk of landslides based on recent precipitation and accumulation of rainwater in the ground.Regarding the format, this is technical but very straightforward. The data format is "General Representation of fields In Binary" (GRIB), and the specific guidelines are given in detail in the World Meteorological Organization's Manual on Codes. The JMA data that we will be examining here follows the FM92 GRIB edition 2, available at the WMO website in English (https://www.wmo.int/pages/prog/www/WMOCodes.html), and in Japanese from the JMA website (http://www.jma.go.jp/jma/kishou/books/tsuhoshiki/tsuhoshiki.html). The Japanese standards are based directly on this FM92 GRIB edition 2, and go by the name 国際気象通報式FM92 GRIB 二進形式格子点資料気象通報式(第2版).There is really no choice but to go through the binary files one byte at a time to check everything out once; then everything can be easily automated. Useful resources when prying through these GRIB2 files are as follows:- JMA Technical Report 374 (配信資料に関する技術情報（気象編）第374号). Available at: http://www.data.jma.go.jp/add/suishin/jyouhou/pdf/374.pdf- Description of the run-length encoding used for compression (ランレングス符号化法の解説). Available at: http://www.jmbsc.or.jp/jp/online/c-onlineGsd.htmlWith these two resources in hand, we open the connection with an example file.<jupyter_code>import os

toread = "demo_data/20130701/Z__C_RJTD_20130701000000_MET_INF_Jdosha_Ggis5km_ANAL_grib2.bin"

f_bin = open(toread, mode="rb")
print(f_bin)<jupyter_output><_io.BufferedReader name='demo_data/20130701/Z__C_RJTD_20130701000000_MET_INF_Jdosha_Ggis5km_ANAL_grib2.bin'><jupyter_text>The files are broken up into sections (節). We shall take one at a time.  Section 0: 指示節This section spans the first sixteen octets.<jupyter_code>f_bin.seek(0)
b = f_bin.read(4) # Should be "GRIB".
print("bytes:", b)
s0_grib = str(b, 'ascii')
print(s0_grib)
b = f_bin.read(2) # Should be missing.
print("bytes:", b)
s0_missing = b
print(s0_missing)<jupyter_output>bytes: b'\xff\xff'
b'\xff\xff'<jupyter_text>As regards the seventh byte, they refer to 符号表0.0, saying that it represents the meteorology field of interest.<jupyter_code>b = f_bin.read(1) # Should be 0.
print("bytes:", b)
s0_field = int.from_bytes(b, byteorder="big", signed=False)
print(s0_field)
b = f_bin.read(1) # Should be 2.
print("bytes:", b)
s0_version = int.from_bytes(b, byteorder="big", signed=False)
print(s0_version)<jupyter_output>bytes: b'\x02'
2<jupyter_text>The final handful of bytes give the *length of the GRIB message* (GRIB報全体の長さ). This is apparently variable from file to file.<jupyter_code>b = f_bin.read(8) # Variable.
print("bytes:", b)
s0_griblen = int.from_bytes(b, byteorder="big", signed=False)
print(s0_griblen)
s0_seclen = 16 # hard-code for convenience.<jupyter_output><empty_output><jupyter_text>___  Section 1: 識別節This section spans a total of twenty-one bytes.<jupyter_code>f_bin.seek(s0_seclen)
b = f_bin.read(4) # Should be 21.
print("bytes:", b)
s1_seclen = int.from_bytes(b, byteorder="big", signed=False)
print(s1_seclen)
b = f_bin.read(1) # Should be 1.
print("bytes:", b)
s1_secnum = int.from_bytes(b, byteorder="big", signed=False)
print(s1_secnum)<jupyter_output>bytes: b'\x01'
1<jupyter_text>The next pair of octects are the "location where the data was created", I assume. They call it “作成中枢の識別”. This is identified via 共通符号表C-1. In any case, the value is supposed to be 34, which apparently corresponds to Tokyo.<jupyter_code>b = f_bin.read(2) # Should be 34.
print("bytes:", b)
s1_center = int.from_bytes(b, byteorder="big", signed=False)
print(s1_center)
b = f_bin.read(2) # Should be 0.
print("bytes:", b)
s1_subcenter = int.from_bytes(b, byteorder="big", signed=False)
print(s1_subcenter)<jupyter_output>bytes: b'\x00\x00'
0<jupyter_text>The next three bytes are rather technical. First is the "version number of GRIB master table" (GRIBマスター表バージョン番号), second is "version number of GRIB region table" (GRIB地域表バージョン番号), and finally is the "meaning of the referred time (参照時刻の意味). Respectively, the meaning of these codes can be found in 符号表1.0, 1.1, and 1.2. Their expected values are 2, 1, and 0.<jupyter_code>b = f_bin.read(1) # Should be 2.
print("bytes:", b)
s1_vermaster = int.from_bytes(b, byteorder="big", signed=False)
print(s1_vermaster)
b = f_bin.read(1) # Should be 1.
print("bytes:", b)
s1_verregion = int.from_bytes(b, byteorder="big", signed=False)
print(s1_verregion)
b = f_bin.read(1) # Should be 0.
print("bytes:", b)
s1_timemeaning = int.from_bytes(b, byteorder="big", signed=False)
print(s1_timemeaning)<jupyter_output>bytes: b'\x00'
0<jupyter_text>Next is some rather important metadata, namely the date/time, which has the meaning specified in the previous byte, which is apparently “解析”. The bytes are as follows:- 2 bytes: year- 1 byte: month- 1 byte: day- 1 byte: hour- 1 byte: minute- 1 byte: secondThe interpretation of these values is specified in 「※１ 時刻表現」. Let's first read them out and take a look.<jupyter_code>b = f_bin.read(2) # Should be the year.
print("bytes:", b)
s1_year = int.from_bytes(b, byteorder="big", signed=False)
print(s1_year)
b = f_bin.read(1) # Should be the month.
print("bytes:", b)
s1_month = int.from_bytes(b, byteorder="big", signed=False)
print(s1_month)
b = f_bin.read(1) # Should be the day.
print("bytes:", b)
s1_day = int.from_bytes(b, byteorder="big", signed=False)
print(s1_day)
b = f_bin.read(1) # Should be the hour.
print("bytes:", b)
s1_hour = int.from_bytes(b, byteorder="big", signed=False)
print(s1_hour)
b = f_bin.read(1) # Should be the minute.
print("bytes:", b)
s1_minute = int.from_bytes(b, byteorder="big", signed=False)
print(s1_minute)
b = f_bin.read(1) # Should be the second.
print("bytes:", b)
s1_second = int.from_bytes(b, byteorder="big", signed=False)
print(s1_second)<jupyter_output>bytes: b'\x00'
0<jupyter_text>Great. Very simple. This looks to match up directly with the file name (i.e., we could just parse the file name to get the same metadata). The next byte is the document's status, either a "product" (現業プロダクト) or "test product" (現業的試験プロジェクト). Details are in 符号表1.3.<jupyter_code>b = f_bin.read(1) # Should be 0 or 1.
print("bytes:", b)
s1_status = int.from_bytes(b, byteorder="big", signed=False)
print(s1_status)<jupyter_output>bytes: b'\x00'
0<jupyter_text>Finally, the document type, detailed in 符号表1.4, is supposedly 2, representing “解析及び予報プロダクト”.<jupyter_code>b = f_bin.read(1) # Should be 2.
print("bytes:", b)
s1_type = int.from_bytes(b, byteorder="big", signed=False)
print(s1_type)<jupyter_output>bytes: b'\x02'
2<jupyter_text>___  Section 2: 地域使用節This section is apparently not used. ___  Section 3: 格子系定義節This section is for lattice definitions, and spans seventy-two bytes.<jupyter_code>f_bin.seek((s0_seclen+s1_seclen))
b = f_bin.read(4) # Should be 72.
print("bytes:", b)
s3_seclen = int.from_bytes(b, byteorder="big", signed=False)
print(s3_seclen)
b = f_bin.read(1) # Should be 3.
print("bytes:", b)
s3_secnum = int.from_bytes(b, byteorder="big", signed=False)
print(s3_secnum)<jupyter_output>bytes: b'\x03'
3<jupyter_text>Reference for the lattice definition, detailed in 符号表3.0. Should be zero, with a meaning of “符号表3.1参照”, rather humorously.<jupyter_code>b = f_bin.read(1) # Should be 0.
print("bytes:", b)
s3_defref = int.from_bytes(b, byteorder="big", signed=False)
print(s3_defref)<jupyter_output>bytes: b'\x00'
0<jupyter_text>Number of points on the lattice.<jupyter_code>b = f_bin.read(4) # Should be 286720.
print("bytes:", b)
s3_numpoints = int.from_bytes(b, byteorder="big", signed=False)
print(s3_numpoints)<jupyter_output>bytes: b'\x00\x04`\x00'
286720<jupyter_text>Number of bytes used for a list that defines the number of points on the lattice (convoluted...), followed by the description of this list. Both are just zero here.<jupyter_code>b = f_bin.read(1) # Should be 0.
print("bytes:", b)
s3_numpointsMem = int.from_bytes(b, byteorder="big", signed=False)
print(s3_numpointsMem)
b = f_bin.read(1) # Should be 0.
print("bytes:", b)
s3_numpointsExp = int.from_bytes(b, byteorder="big", signed=False)
print(s3_numpointsExp)<jupyter_output>bytes: b'\x00'
0<jupyter_text>The next two bytes give the template number identifying the lattice definition type, details in 符号表3.1. It should be zero, meaning “緯度・経度格子”.<jupyter_code>b = f_bin.read(2) # Should be 0.
print("bytes:", b)
s3_template = int.from_bytes(b, byteorder="big", signed=False)
print(s3_template)<jupyter_output>bytes: b'\x00\x00'
0<jupyter_text>Next is a key value: the identifier of the shape of the earth, details in 符号表3.2. This should be four, representing “GRS80回転楕円体”, recalling GRS80 is *Geodetic Reference System 1980*, a reference ellipsoid modelling the earth.<jupyter_code>b = f_bin.read(1) # Should be 4.
print("bytes:", b)
s3_earthshape = int.from_bytes(b, byteorder="big", signed=False)
print(s3_earthshape)<jupyter_output>bytes: b'\x04'
4<jupyter_text>The next five are all missing.<jupyter_code>b = f_bin.read(5) # Should be missing, likely max values.
print("bytes:", b)<jupyter_output>bytes: b'\xff\xff\xff\xff\xff'<jupyter_text>The next ten bytes are broken up as follows.- 1 byte: major axis unit factor (尺度因子)- 4 bytes: major axis length in given units (長軸の尺度付きの長さ)- 1 byte: minor axis unit factor (尺度因子)- 4 bytes: minor axis length in given units (短軸の尺度付きの長さ)<jupyter_code>b = f_bin.read(1) # Should be 1.
print("bytes:", b)
s3_factormaj = int.from_bytes(b, byteorder="big", signed=False)
print(s3_factormaj)
b = f_bin.read(4) # Should be 63781370.
print("bytes:", b)
s3_lenmaj = int.from_bytes(b, byteorder="big", signed=False)
print(s3_lenmaj)
b = f_bin.read(1) # Should be 1.
print("bytes:", b)
s3_factormin = int.from_bytes(b, byteorder="big", signed=False)
print(s3_factormin)
b = f_bin.read(4) # Should be 63567523.
print("bytes:", b)
s3_lenmin = int.from_bytes(b, byteorder="big", signed=False)
print(s3_lenmin)<jupyter_output>bytes: b'\x03\xc9\xf6\xa3'
63567523<jupyter_text>Next, we have the absolute number of points on the lattice, following along the latitude and longitude lines. In Japanese, the former is “緯線に沿った格子点数”, the latter is “経線に沿った格子点数”.<jupyter_code>b = f_bin.read(4) # Should be 512.
print("bytes:", b)
s3_numlon = int.from_bytes(b, byteorder="big", signed=False)
print(s3_numlon)
b = f_bin.read(4) # Should be 560.
print("bytes:", b)
s3_numlat = int.from_bytes(b, byteorder="big", signed=False)
print(s3_numlat)<jupyter_output>bytes: b'\x00\x00\x020'
560<jupyter_text>The interpretation as a grid is easy: we have a grid which is 512 cells wide, and 560 cells tall. As well, the number of cells computed here exactly equals the total number of points on the lattice we extracted earlier.<jupyter_code>print(s3_numpoints)
print((s3_numlat*s3_numlon))<jupyter_output>286720
286720<jupyter_text>Next is some "base angle", should be zero.<jupyter_code>b = f_bin.read(4) # Should be 0.
print("bytes:", b)
s3_baseangle = int.from_bytes(b, byteorder="big", signed=False)
print(s3_baseangle)<jupyter_output>bytes: b'\x00\x00\x00\x00'
0<jupyter_text>Next should be missing.<jupyter_code>b = f_bin.read(4) # Should be 0.
print("bytes:", b)<jupyter_output>bytes: b'\xff\xff\xff\xff'<jupyter_text>Now are some numerical values we need to be careful with. The next eight bytes give the latitude (緯度) and longitude (経度) of the "first" point on the lattice. The units are in what they call “10-6度単位”, which is scientific notation for units re-scaled by $10^{-6}$. Let's see the values first, and then confirm we are handling what we are expecting.<jupyter_code>b = f_bin.read(4) # Should be 47975000.
print("bytes:", b)
s3_firstlat = int.from_bytes(b, byteorder="big", signed=False)
print(s3_firstlat)
b = f_bin.read(4) # Should be 118031250.
print("bytes:", b)
s3_firstlon = int.from_bytes(b, byteorder="big", signed=False)
print(s3_firstlon)<jupyter_output>bytes: b'\x07\t\x03\x92'
118031250<jupyter_text>Note that if we multiply these values by $10^{-6}$, we have:<jupyter_code>print((s3_firstlat/10**6))
print((s3_firstlon/10**6))<jupyter_output>47.975
118.03125<jupyter_text>This certainly look a lot more like degrees, and of course that is precisely what they are. The detailed calculations given in the JMA materials are as follows:- First point's latitude: 48N - 3.0/60/2.- First point's longitude: 118E + 3.75/60/2.Of course, these calculations are easily confirmed:<jupyter_code>print((48 - (3/(2*60))))
print((118 + (3.75/(2*60))))<jupyter_output>47.975
118.03125<jupyter_text>Next is a slightly complicated pair of terms: 分解能及び成分フラグ, with details in フラグ表3.3. The value should be 0x30, though it comes up as zero. In any case, the values after this guy are fine, so it seems there is no big issue.<jupyter_code>b = f_bin.read(1) # Should be 0x30.
print("bytes:", b)<jupyter_output>bytes: b'0'<jupyter_text>Next, just as before, we have eight bytes dedicated to the latitude and longitude of the final lattice point/<jupyter_code>b = f_bin.read(4) # Should be 20025000.
print("bytes:", b)
s3_lastlat = int.from_bytes(b, byteorder="big", signed=False)
print(s3_lastlat)
b = f_bin.read(4) # Should be 149968750.
print("bytes:", b)
s3_lastlon = int.from_bytes(b, byteorder="big", signed=False)
print(s3_lastlon)<jupyter_output>bytes: b'\x08\xf0Wn'
149968750<jupyter_text>We can do analogous confirmations as before, although noting that the +/- signs in the shift are now opposite what they were for the first point.- Last point's latitude: 20N + 3.0/60/2.- Last point's longitude: 150E - 3.75/60/2.Confirming everything:<jupyter_code>print((s3_lastlat/10**6))
print((s3_lastlon/10**6))
print((20 + (3/(2*60))))
print((150 - (3.75/(2*60))))<jupyter_output>20.025
149.96875<jupyter_text>Finally, we have the step size in the $i$-direction and $j$-direction respectively (called the $i$方向の増分, $j$方向の増分), defined:- $i$-direction: 3.75/60- $j$-direction: 3.0/60The units, again, are $10^{-6}$ as before.<jupyter_code>b = f_bin.read(4) # Should be 62500.
print("bytes:", b)
s3_isize = int.from_bytes(b, byteorder="big", signed=False)
print(s3_isize)
b = f_bin.read(4) # Should be 50000.
print("bytes:", b)
s3_jsize = int.from_bytes(b, byteorder="big", signed=False)
print(s3_jsize)
print((s3_isize/10**6))
print((s3_jsize/10**6))
print(3.75/60)
print(3.0/60)<jupyter_output>0.0625
0.05<jupyter_text>As we would expect, the degrees *increase uniformly* (by the $i$ and $j$ sizes) over the lat/long directions on our graph.<jupyter_code>tmp_val = s3_firstlon
print(tmp_val)
for i in range(s3_numlon-1):
    tmp_val += s3_isize
print(tmp_val)
print(s3_lastlon)
tmp_val = s3_firstlat
print(tmp_val)
for j in range(s3_numlat-1):
    tmp_val -= s3_jsize
print(tmp_val)
print(s3_lastlat)<jupyter_output>47975000
20025000
20025000<jupyter_text>Finally, we have the “走査モード”, namely a scanning mode, details in フラグ表3.4, and it should have a value of 0x00, namely zero.<jupyter_code>b = f_bin.read(1) # Should be 0x00.
print("bytes:", b)<jupyter_output>bytes: b'\x00'<jupyter_text>___  Section 4: プロダクト定義節This section is for so-called "product" definitions, and spans forty-two bytes.<jupyter_code>f_bin.seek((s0_seclen+s1_seclen+0+s3_seclen))
b = f_bin.read(4) # Should be 42.
print("bytes:", b)
s4_seclen = int.from_bytes(b, byteorder="big", signed=False)
print(s4_seclen)
b = f_bin.read(1) # Should be 4.
print("bytes:", b)
s4_secnum = int.from_bytes(b, byteorder="big", signed=False)
print(s4_secnum)<jupyter_output>bytes: b'\x04'
4<jupyter_text>The number of coordinate values immediately after the template (テンプレート直後の座標値の数). Should be zero.<jupyter_code>b = f_bin.read(2) # Should be 0.
print("bytes:", b)
s4_coordnum = int.from_bytes(b, byteorder="big", signed=False)
print(s4_coordnum)<jupyter_output>bytes: b'\x00\x00'
0<jupyter_text>Product definition template number. Should be 50000, apparently with the meaning that the "product" here is based on processing other existing "products". Details in 符号表4.0.<jupyter_code>b = f_bin.read(2) # Should be 50000.
print("bytes:", b)
s4_tempnum = int.from_bytes(b, byteorder="big", signed=False)
print(s4_tempnum)<jupyter_output>bytes: b'\xc3P'
50000<jupyter_text>Parameter category. Should be one, which means humidity (湿度), details in 符号表4.1.<jupyter_code>b = f_bin.read(1) # Should be 1.
print("bytes:", b)
s4_paracat = int.from_bytes(b, byteorder="big", signed=False)
print(s4_paracat)<jupyter_output>bytes: b'\x01'
1<jupyter_text>Parameter number. Should be 208, which means landslide risk warnings (土砂災害警戒判定値), details in 符号表4.2.<jupyter_code>b = f_bin.read(1) # Should be 208.
print("bytes:", b)
s4_paranum = int.from_bytes(b, byteorder="big", signed=False)
print(s4_paranum)<jupyter_output>bytes: b'\xd0'
208<jupyter_text>Type of processing, details in 符号表4.3. Should be zero, meaning analysis/forecasting (解析及び予報).<jupyter_code>b = f_bin.read(1) # Should be 0.
print("bytes:", b)
s4_proctype = int.from_bytes(b, byteorder="big", signed=False)
print(s4_proctype)<jupyter_output>bytes: b'\x00'
0<jupyter_text>Background processing method, details in 符号表JMA4.1, should be 160, meaning the specialized routine for landslide risk warnings (土砂災害警戒情報ルーチン).<jupyter_code>b = f_bin.read(1) # Should be 160.
print("bytes:", b)
s4_backproc = int.from_bytes(b, byteorder="big", signed=False)
print(s4_backproc)<jupyter_output>bytes: b'\xa0'
160<jupyter_text>Next one is missing, the analysis method identifier.<jupyter_code>b = f_bin.read(1) # Should be missing.
print("bytes:", b)<jupyter_output>bytes: b'\xff'<jupyter_text>Cut-off time measured from the reference time, in hours and minutes.<jupyter_code>b = f_bin.read(2) # Should be 0.
print("bytes:", b)
s4_cutoffHour = int.from_bytes(b, byteorder="big", signed=False)
print(s4_cutoffHour)
b = f_bin.read(1) # Should be 10.
print("bytes:", b)
s4_cutoffMin = int.from_bytes(b, byteorder="big", signed=False)
print(s4_cutoffMin)<jupyter_output>bytes: b'\n'
10<jupyter_text>Type of units used for time period, details in 符号表4.4. Should be zero, meaning minutes.<jupyter_code>b = f_bin.read(1) # Should be 0.
print("bytes:", b)
s4_timeunit = int.from_bytes(b, byteorder="big", signed=False)
print(s4_timeunit)<jupyter_output>bytes: b'\x00'
0<jupyter_text>Forecast time (予報時間). This is supposed to be zero.<jupyter_code>b = f_bin.read(4) # Should be 0.
print("bytes:", b)
s4_fctime = int.from_bytes(b, byteorder="big", signed=False)
print(s4_fctime)<jupyter_output>bytes: b'\x00\x00\x00\x00'
0<jupyter_text>Type of "first fixed surface" (第一固定面の種類), with details in 符号表4.5. Should be one, meaning either land or water (地面または水面).<jupyter_code>b = f_bin.read(1) # Should be 1.
print("bytes:", b)
s4_surftype = int.from_bytes(b, byteorder="big", signed=False)
print(s4_surftype)<jupyter_output>bytes: b'\x01'
1<jupyter_text>The next eleven bytes are all missing.<jupyter_code>b = f_bin.read(1) # Should be missing.
print("bytes:", b)
b = f_bin.read(4) # Should be missing.
print("bytes:", b)
b = f_bin.read(1) # Should be missing.
print("bytes:", b)
b = f_bin.read(1) # Should be missing.
print("bytes:", b)
b = f_bin.read(4) # Should be missing.
print("bytes:", b)<jupyter_output>bytes: b'\xff\xff\xff\xff'<jupyter_text>Name of the first set of related materials used in making this data set, details in 符号表JMA4.5, should be 2, meaning 土壌雨量指数の予測値. Otherwise, if the values are missing, it will be 255.<jupyter_code>b = f_bin.read(1) # Should be 2 or 255.
print("bytes:", b)
s4_refsoil = int.from_bytes(b, byteorder="big", signed=False)
print(s4_refsoil)<jupyter_output>bytes: b'\x02'
2<jupyter_text>Difference between the *analysis* time and the reference time in the *current* referenced materials, in hours and minutes. By "current", we mean the materials identified in the previous octet. After this, there will be another set of reference materials, and of course in principle this could go on indefinitely, but here just two reference materials in total.<jupyter_code>b = f_bin.read(2) # Should be variable.
print("bytes:", b)
s4_refsoilDiffHour = int.from_bytes(b, byteorder="big", signed=False)
print(s4_refsoilDiffHour)
b = f_bin.read(1) # Should be variable.
print("bytes:", b)
s4_refsoilDiffMin = int.from_bytes(b, byteorder="big", signed=False)
print(s4_refsoilDiffMin)<jupyter_output>bytes: b'&'
38<jupyter_text>Name of the second set of related materials used in making this data set, details in 符号表JMA4.5, should be 4, meaning 1時間降水量の予測値. Otherwise, if the values are missing, it will be 255.<jupyter_code>b = f_bin.read(1) # Should be 4 or 255.
print("bytes:", b)
s4_refrain = int.from_bytes(b, byteorder="big", signed=False)
print(s4_refrain)<jupyter_output>bytes: b'\x04'
4<jupyter_text>Once again, time differences, now for the second set of reference materials.<jupyter_code>b = f_bin.read(2) # Should be variable.
print("bytes:", b)
s4_refrainDiffHour = int.from_bytes(b, byteorder="big", signed=False)
print(s4_refrainDiffHour)
b = f_bin.read(1) # Should be variable.
print("bytes:", b)
s4_refrainDiffMin = int.from_bytes(b, byteorder="big", signed=False)
print(s4_refrainDiffMin)<jupyter_output>bytes: b'\x00'
0<jupyter_text>For reference, from the tenth byte through to the final byte of this section above, they say that "template 4.50000" is used. ___  Section 5: 資料表現節This section has "representations" of the main "materials". This is not the actual data, but rather the means of interpreting the data values which appear in the seventh section, a bit later.<jupyter_code>f_bin.seek((s0_seclen+s1_seclen+0+s3_seclen+s4_seclen))
b = f_bin.read(4) # Should be variable.
print("bytes:", b)
s5_seclen = int.from_bytes(b, byteorder="big", signed=False)
print(s5_seclen)
b = f_bin.read(1) # Should be 5.
print("bytes:", b)
s5_secnum = int.from_bytes(b, byteorder="big", signed=False)
print(s5_secnum)<jupyter_output>bytes: b'\x05'
5<jupyter_text>Total number of data instances, should be equal to the number of lattice points seen earlier.<jupyter_code>b = f_bin.read(4) # Should be 286720.
print("bytes:", b)
s5_numdata = int.from_bytes(b, byteorder="big", signed=False)
print(s5_numdata)<jupyter_output>bytes: b'\x00\x04`\x00'
286720<jupyter_text>Next is the template number of these materials, details in 符号表5.0. Should be 200, with meaning of 格子点資料 with run-length compression (ランレングス圧縮).<jupyter_code>b = f_bin.read(2) # Should be 200.
print("bytes:", b)
s5_tempnum = int.from_bytes(b, byteorder="big", signed=False)
print(s5_tempnum)<jupyter_output>bytes: b'\x00\xc8'
200<jupyter_text>Next is the number of bits per individual datum.<jupyter_code>b = f_bin.read(1) # Should be 8.
print("bytes:", b)
s5_databits = int.from_bytes(b, byteorder="big", signed=False)
print(s5_databits)<jupyter_output>bytes: b'\x08'
8<jupyter_text>Maximum "level" used in compression this time, a potentially variable value. Defined $V$ for later algebraic formulations.<jupyter_code>b = f_bin.read(2) # Should be variable.
print("bytes:", b)
s5_V = int.from_bytes(b, byteorder="big", signed=False)
print(s5_V)<jupyter_output>bytes: b'\x00\x04'
4<jupyter_text>Maximum level, should be ten, and defined as $M$ for later use. Also, should have $V \leq M$.<jupyter_code>b = f_bin.read(2) # Should be 10.
print("bytes:", b)
s5_M = int.from_bytes(b, byteorder="big", signed=False)
print(s5_M)<jupyter_output>bytes: b'\x00\n'
10<jupyter_text>Representative data value unit factor, should be zero.<jupyter_code>b = f_bin.read(1) # Should be 0.
print("bytes:", b)
s5_datafactor = int.from_bytes(b, byteorder="big", signed=False)
print(s5_datafactor)<jupyter_output>bytes: b'\x00'
0<jupyter_text>Now, ranging over $m=1,2,\ldots,M$, we look at the data values that correspond to level $m$ (レベル$m$に対応するデータ代表値). Each of these uses two bytes. The interpretation is as follows (via ※2 in documentation):- $m=0$. $R(m)$ is null. This is missing value.- $m=1$. $R(m)=-2$. This is ocean or non-JPN lattice point.- $m=2$. $R(m)=-1$. Lattice point on Japanese land, but not applicable for warnings.- $m=3$. $R(m)=0$. This is warning level 0.- $m=4$. $R(m)=1$. This is warning level 1.- $m=5$. $R(m)=2$. This is warning level 2.- $m=6$. $R(m)=3$. This is warning level 3.- $m=7$. $R(m)=4$. This is warning level 4.- $m=8$. $R(m)=5$. This is TBD.- $m=9$. $R(m)=6$. This is TBD.- $m=10$. $R(m)=7$. This is TBD.Let's take a look at the first guys, one byte at a time.<jupyter_code>b = f_bin.read(1)
print("bytes:", b)
b = f_bin.read(1)
print("bytes:", b)
b = f_bin.read(1)
print("bytes:", b)
b = f_bin.read(1)
print("bytes:", b)<jupyter_output>bytes: b'\x80'
bytes: b'\x02'
bytes: b'\x80'
bytes: b'\x01'<jupyter_text>Clearly, they're just using the byte value (in hex) of 0x80 to be a minus sign, as the second byte in each pair is precisely as we expect (2 and 1). Checking out the rest is easy.<jupyter_code>for m in range(3,(s5_M+1)):
    b = f_bin.read(2)
    print("bytes:", b)
    tmp_level = int.from_bytes(b, byteorder="big", signed=True)
    print("m =", m, "; R(m) =", tmp_level)<jupyter_output>bytes: b'\x00\x00'
m = 3 ; R(m) = 0
bytes: b'\x00\x01'
m = 4 ; R(m) = 1
bytes: b'\x00\x02'
m = 5 ; R(m) = 2
bytes: b'\x00\x03'
m = 6 ; R(m) = 3
bytes: b'\x00\x04'
m = 7 ; R(m) = 4
bytes: b'\x00\x05'
m = 8 ; R(m) = 5
bytes: b'\x00\x06'
m = 9 ; R(m) = 6
bytes: b'\x00\x07'
m = 10 ; R(m) = 7<jupyter_text>___  Section 6: ビットマップ節This section is basically empty, spanning just six bytes.<jupyter_code>f_bin.seek((s0_seclen+s1_seclen+0+s3_seclen+s4_seclen+s5_seclen))
b = f_bin.read(4) # Should be 6.
print("bytes:", b)
s6_seclen = int.from_bytes(b, byteorder="big", signed=False)
print(s6_seclen)
b = f_bin.read(1) # Should be 6.
print("bytes:", b)
s6_secnum = int.from_bytes(b, byteorder="big", signed=False)
print(s6_secnum)<jupyter_output>bytes: b'\x06'
6<jupyter_text>Bitmap directive, should be 255, meaning that a bitmap isn't used here.<jupyter_code>b = f_bin.read(1) # Should be 255.
print("bytes:", b)
s6_bitmap = int.from_bytes(b, byteorder="big", signed=False)
print(s6_bitmap)<jupyter_output>bytes: b'\xff'
255<jupyter_text>___  Section 7: 資料節Finally, the data section.<jupyter_code>f_bin.seek((s0_seclen+s1_seclen+0+s3_seclen+s4_seclen+s5_seclen+s6_seclen))
b = f_bin.read(4) # Should be variable.
print("bytes:", b)
s7_seclen = int.from_bytes(b, byteorder="big", signed=False)
print(s7_seclen)
b = f_bin.read(1) # Should be 7.
print("bytes:", b)
s7_secnum = int.from_bytes(b, byteorder="big", signed=False)
print(s7_secnum)<jupyter_output>bytes: b'\x07'
7<jupyter_text>Aside: JMA run-length compression algorithmNow comes a long series of data *values*, although these are "run-length compressed octets". The decompression is extremely straightforward once we know how the compression is done. Fortunately, this is explained in the second reference material linked at the start of this notebook, and we spell it out here for clarity and posterity.Run-length encoding is a very simple (lossless) compression technique. Given a sequence of scalar values $(d_{1},d_{2},\ldots,d_{k})$, if there are long stretches where values repeat, i.e. $d_{i}=d_{i+1}=\cdots=d_{i+m}$, then the basic idea is to just store (1) the value (here, $d_{i}$), and (2) the "run length" (here, $m+1$). With this information, we can naturally re-construct the original sequence.The JMA run-length encoding for domestic meteorological data stored in binary format is just one special case of run length encoding, and we describe it in detail now.First, two key terms to get us started.- `MAXV` (here $V$): all data values $d$ are assumed to be integers $0 \leq d \leq V$.- `NBIT` (here $N$): the number of bits that makes up a single data point in the compressed data.Now let's talk about the content of these data points. During compression, data are organized into "sets" (contiguous sub-sequences). There are precisely two types of sets: - Sets with just one data point $(c)$, holding *value* information. - Sets with two or more data points $(c,r_{0},r_{1},\ldots)$, holding *value* and *run length* information.  Each of these data points is provided in the compressed data set using a base-2 (binary) representation, composed of `NBIT` binary digits. In both cases, $c$ represents the actual data value of interest, but the former case indicates that the value appeared just *once*, whereas in the latter case, the $r_{0},r_{1},\ldots$ provide us all the information needed to determine the run length (by construction greater than $1$, details given below).Given compressed data $(c_{1},\ldots,c_{m})$, determining where one "set" ends and another begins is simple, using the following rules:- If $0 \leq c_{i} \leq V$, then $c_{i}$ is a *value*, not run length information.- Else, we have $c_{i} > V$ and $c_{i}$ is run length information for a previous value.To completely decode this data, we need not just $c$, but the actual run length, let us denote this by $R$. In the case of the trivial set $(c)$, we know $R=1$. In general, however, given a set $(c,r_{0},r_{1},\ldots)$, we need to deduce $R$ from the $r_{0},r_{1},\ldots$ values. We describe how to do this next.Since our compressed data has only `NBIT` bits for each datum, naturally each datum can represent at most $2^{N}$ different values. On the other hand, run lengths are at least $1$, but can assuredly be much larger than $2^{N}$. This means we will need more than one datum. Furthermore, for very long runs, if we were to simply sum the $r_{0},r_{1},\ldots$ to get the run length, this would require a lot of data, again very inefficient.We have the following requirements:- We would like to represent run length $R$ using as little data as possible.- Each $r_{i}$ must satisfy $r_{i}>V$ in order to distinguish run length information from value information.The technique employed is to use a base-$L$ (written `LNGU` in JMA materials and our code) representation of run length, where $L=2^{N}-V-1$. Write $a_{p-1}a_{p-2} \cdots a_{0}$ to denote a $p$-digit base-$L$ representation of $R$, namely\begin{align*}R = a_{0}L^{0} + \cdots + a_{p-1}L^{p-1}\end{align*}As each digit $a_{j}$ ranges between $0$ and $L-1 = 2^{N}-V-2$, we cannot simply use these $a_{j}$ values as our $r_{i}$ values, since they cannot be distinguished from data values in general. This problem is easily solved: pad the values by adding $V+1$, such that they range between $V+1$ and $(L-1)+(V+1)=2^{N}-1$. This upper bound is precisely the maximum value that we can represent using `NBIT` bits, and the lower bound is greater than `MAXV`, as required.Written out explicitly, we have $r_{i} = a_{i} + V + 1$. This is all the information we need to decode the run length.In practice, we start with the $c$ and $r_{i}$ values, so we have to work backwards:0. Given $V$, and $N$ (parameters of the compression procedure), we can compute $L$.0. If given a set $(c)$, decompress as-is with $(c)$.0. Else, given a set $(c,r_{0},\ldots,r_{p-1})$, decompress to $(c,\ldots,c)$ repeated $R > 1$ times, with $R$ computed as follows:   0. Compute $a_{i} = r_{i}-V-1$ for $i=0,\ldots,p-1$.   0. Given $L$ and $(a_{0},\ldots,a_{p-1})$, compute $R$ as $a_{p-1} \cdots a_{0}$ in base-$L$.__On data orientation:__ For 2D lattice data stored in a one-dimensional array (a long vector), the assumption is that it has been read going from *west to east*, and *north to south*.  Decompression of lattice dataNow that we've covered the background of the compression algorithm used to create the data we have accessto, moving forward we read the key parameters, then deploy the simple decompression algorithm described above.<jupyter_code>bytes_left = s7_seclen - 5 # number of bytes left in this section.
# Before doing anything, read out the compressed data sequence.
data_compressed = []
bpd = s5_databits // 8 # bytes per data point.
for i in range(bytes_left):
    b = f_bin.read(bpd)
    #print("bytes:", b)
    b_int = int.from_bytes(b, byteorder="big", signed=False)
    data_compressed += [b_int]
# Key constants for decompression.
NBIT = s5_databits
MAXV = s5_V
LNGU = 2**NBIT - 1 - MAXV<jupyter_output><empty_output><jupyter_text>The `NBIT`, `MAXV`, and `LNGU` values for this particular data set are as follows.<jupyter_code>print(NBIT, MAXV, LNGU)
import numpy as np
# Next, decompress this data sequence.
data_decompressed = []
dictpair = {"data": [], "runlen": []}
for i in range(len(data_compressed)):
    
    tocheck = data_compressed[i]
    print("tocheck =", tocheck, "; dict =", dictpair)
    
    if len(dictpair["data"]) == 0:
        dictpair["data"] = [tocheck]
    else:
        # If data slot is populated, then we're either dealing
        # with runlength, or moving to a new singleton.
        
        if tocheck > MAXV:
            # In this case, it is runlength data to be added.
            dictpair["runlen"] += [tocheck]
        else:
            # In this case, it is data, signifying a new set.
            # Thus, we must enter the "processing" routine.
            
            numrl = len(dictpair["runlen"])
            if numrl == 0:
                # If no runlength info, means it only appeared once.
                data_decompressed += dictpair["data"]
                #print("Adding to decompressed data:", dictpair["data"])
                
            else:
                # If there is runlength info, then must process it.
                rlinfo = np.array(dictpair["runlen"])
                rlexp = np.arange(numrl)
                rlvec = (LNGU**rlexp) * (rlinfo - (1+MAXV))
                rl = np.sum(rlvec) + 1
                data_decompressed += rl * dictpair["data"]
                #print("Adding to decompressed data:", rl * dictpair["data"])
            
            # Reset, critically noting that we have a new data point to store.
            dictpair["data"] = [tocheck]
            dictpair["runlen"] = []

# Finally, note that there will ALWAYS be something left over.
if tocheck <= MAXV:
    # If the last guy was small, then just add it.
    data_decompressed += dictpair["data"]
    #print("Adding to decompressed data:", dictpair["data"])
else:
    # If the last guy was large, we have non-trivial run-length.
    numrl = len(dictpair["runlen"])
    rlinfo = np.array(dictpair["runlen"])
    rlexp = np.arange(numrl)
    rlvec = (LNGU**rlexp) * (rlinfo - (1+MAXV))
    rl = np.sum(rlvec) + 1
    data_decompressed += rl * dictpair["data"]
    #print("Adding to decompressed data:", rl * dictpair["data"])
print(len(data_decompressed))
print(np.unique(np.array(data_decompressed)))
#print(data_decompressed)
print(s3_numpoints)
#print(512*560)
#print(s3_numlat)
#print(s3_numlon)
s3_numpoints-len(data_decompressed)<jupyter_output><empty_output><jupyter_text>If this value is zero, we can be essentially certain that we have successfully decoded the compressed data. ___  Section 8: 終端節Basically a trivial final section that lets us know we're finished.<jupyter_code>f_bin.seek((s0_seclen+s1_seclen+0+s3_seclen+s4_seclen+s5_seclen+s6_seclen+s7_seclen))
b = f_bin.read(4) # Should be "7777".
print(b)
s8_grib = str(b, 'ascii')
print(s8_grib)
b = f_bin.read(1) # Should be at the end...
print(b)<jupyter_output>b''
StarcoderdataJupyterStructuredCleanDedup