GRIB学习笔记：样例文件分析 - 位图
目录
本文通过分析一个GRIB2消息介绍使用Section 6位图段的GRIB2数据。
介绍
某些数据可能有缺失值，我们一般使用一个比较大的值代表缺失，如果将该值保存到数据中，会导致存储空间浪费。
Section 6位图段用于处理这种情况，位图段使用1个位表示网格点对应的数据是否保存到数据段中。
- 1：数据保存
- 0：数据不保存
相比于数据段每个数据使用10个bits保存，位图段使用1个bit表示数据缺失能有效减少存储空间。
示例
样例文件来自中国气象局数值预报中心GRAPES全球模式的预报结果，选用2019年11月4日00时次003时效的反照率（Albedo (%) ）。
***** FILE: bitmap.grib2
#============== MESSAGE 1 ( length=600630 ) ==============
1-4 identifier = GRIB
5-6 reserved = MISSING
7 discipline = 0 [Meteorological products (grib2/tables/4/0.0.table) ]
8 editionNumber = 2
9-16 totalLength = 600630
====================== SECTION_1 ( length=21, padding=0 ) ======================
1-4 section1Length = 21
5 numberOfSection = 1
6-7 centre = 38 [Beijing (RSMC) (common/c-11.table) ]
8-9 subCentre = 0
10 tablesVersion = 4 [Version implemented on 7 November 2007 (grib2/tables/1.0.table) ]
11 localTablesVersion = 1 [Unknown code table entry (grib2/tables/4/1.1.table) ]
12 significanceOfReferenceTime = 1 [Start of forecast (grib2/tables/4/1.2.table) ]
13-14 year = 2019
15 month = 11
16 day = 4
17 hour = 0
18 minute = 0
19 second = 0
20 productionStatusOfProcessedData = 0 [Operational products (grib2/tables/4/1.3.table) ]
21 typeOfProcessedData = 1 [Forecast products (grib2/tables/4/1.4.table) ]
====================== SECTION_3 ( length=72, padding=0 ) ======================
1-4 section3Length = 72
5 numberOfSection = 3
6 sourceOfGridDefinition = 0 [Specified in Code table 3.1 (grib2/tables/4/3.0.table) ]
7-10 numberOfDataPoints = 1036800
11 numberOfOctectsForNumberOfPoints = 0
12 interpretationOfNumberOfPoints = 0 [There is no appended list (grib2/tables/4/3.11.table) ]
13-14 gridDefinitionTemplateNumber = 0 [Latitude/longitude (Also called equidistant cylindrical, or Plate Carree) (grib2/tables/4/3.1.table) ]
15 shapeOfTheEarth = 6 [Earth assumed spherical with radius of 6,371,229.0 m (grib2/tables/4/3.2.table) ]
16 scaleFactorOfRadiusOfSphericalEarth = 0
17-20 scaledValueOfRadiusOfSphericalEarth = 0
21 scaleFactorOfEarthMajorAxis = 0
22-25 scaledValueOfEarthMajorAxis = 0
26 scaleFactorOfEarthMinorAxis = 0
27-30 scaledValueOfEarthMinorAxis = 0
31-34 Ni = 1440
35-38 Nj = 720
39-42 basicAngleOfTheInitialProductionDomain = 0
43-46 subdivisionsOfBasicAngle = 0
47-50 latitudeOfFirstGridPoint = 89875000
51-54 longitudeOfFirstGridPoint = 0
55 resolutionAndComponentFlags = 48 [00110000]
56-59 latitudeOfLastGridPoint = -89875000
60-63 longitudeOfLastGridPoint = 359750000
64-67 iDirectionIncrement = 250000
68-71 jDirectionIncrement = 250000
72 scanningMode = 0 [00000000]
====================== SECTION_4 ( length=34, padding=0 ) ======================
1-4 section4Length = 34
5 numberOfSection = 4
6-7 NV = 0
8-9 productDefinitionTemplateNumber = 0 [Analysis or forecast at a horizontal level or in a horizontal layer at a point in time (grib2/tables/4/4.0.table) ]
10 parameterCategory = 19 [Physical atmospheric properties (grib2/tables/4/4.1.0.table) ]
11 parameterNumber = 1 [Albedo (%) (grib2/tables/4/4.2.0.19.table) ]
12 typeOfGeneratingProcess = 2 [Forecast (grib2/tables/4/4.3.table) ]
13 backgroundProcess = 0
14 generatingProcessIdentifier = 15
15-16 hoursAfterDataCutoff = 0
17 minutesAfterDataCutoff = 0
18 indicatorOfUnitOfTimeRange = 1 [Hour (grib2/tables/4/4.4.table) ]
19-22 forecastTime = 3
23 typeOfFirstFixedSurface = 1 [Ground or water surface (grib2/tables/4/4.5.table) ]
24 scaleFactorOfFirstFixedSurface = MISSING
25-28 scaledValueOfFirstFixedSurface = MISSING
29 typeOfSecondFixedSurface = 255 [Missing (grib2/tables/4/4.5.table) ]
30 scaleFactorOfSecondFixedSurface = MISSING
31-34 scaledValueOfSecondFixedSurface = MISSING
====================== SECTION_5 ( length=23, padding=0 ) ======================
1-4 section5Length = 23
5 numberOfSection = 5
6-9 numberOfValues = 500443
10-11 dataRepresentationTemplateNumber = 40 [JPEG2000 Packing (grib2/tables/4/5.0.table) ]
12-15 referenceValue = 0
16-17 binaryScaleFactor = 0
18-19 decimalScaleFactor = 2
20 bitsPerValue = 13
21 typeOfOriginalFieldValues = 0 [Floating point (grib2/tables/4/5.1.table) ]
22 typeOfCompressionUsed = 0 [Lossless (grib2/tables/4/5.40.table) ]
23 targetCompressionRatio = 255
====================== SECTION_6 ( length=129606, padding=0 ) ======================
1-4 section6Length = 129606
5 numberOfSection = 6
6 bitMapIndicator = 0 [A bit map applies to this product and is specified in this Section (grib2/tables/4/6.0.table) ]
7-129606 bitmap = 129600 {
00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00,
00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00,
00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00,
00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00,
00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00,
00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00, 00,
00, 00, 00, 00
... 129500 more values
} # g2bitmap bitmap
====================== SECTION_7 ( length=470854, padding=0 ) ======================
1-4 section7Length = 470854
5 numberOfSection = 7
6-470854 codedValues = (500443,470849) {
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 6.7470000000e+01, 6.7540000000e+01, 6.7490000000e+01,
6.7770000000e+01, 6.7430000000e+01, 6.7830000000e+01, 6.7520000000e+01, 6.7510000000e+01, 6.7660000000e+01, 6.7370000000e+01, 6.7340000000e+01,
6.7480000000e+01, 6.7840000000e+01, 6.7880000000e+01, 6.7950000000e+01, 6.7780000000e+01, 6.7830000000e+01, 6.8030000000e+01, 6.7880000000e+01,
6.8100000000e+01, 6.7840000000e+01, 6.8080000000e+01, 6.8210000000e+01, 6.8120000000e+01, 6.7930000000e+01, 6.8020000000e+01, 6.7980000000e+01,
6.8030000000e+01, 6.8120000000e+01, 6.8310000000e+01, 6.8560000000e+01, 6.8580000000e+01, 6.8910000000e+01, 6.9460000000e+01, 6.9690000000e+01,
6.9490000000e+01, 6.9620000000e+01, 6.9990000000e+01, 7.0100000000e+01, 6.9940000000e+01, 6.9690000000e+01, 6.9950000000e+01, 7.0270000000e+01,
7.0030000000e+01, 6.9880000000e+01, 7.0100000000e+01, 7.0410000000e+01, 7.0100000000e+01, 6.9570000000e+01, 6.9340000000e+01, 0.0000000000e+00,
0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 6.8140000000e+01, 6.7800000000e+01, 6.7810000000e+01, 6.8310000000e+01,
6.8290000000e+01, 6.7820000000e+01, 6.7490000000e+01, 6.7660000000e+01, 6.7420000000e+01, 6.7890000000e+01, 6.7000000000e+01, 6.6680000000e+01,
6.6530000000e+01, 6.6570000000e+01, 6.6720000000e+01, 6.6640000000e+01, 6.6740000000e+01, 6.6920000000e+01, 6.7200000000e+01, 6.7360000000e+01,
6.7330000000e+01, 6.7080000000e+01, 6.7630000000e+01, 6.7650000000e+01, 6.7560000000e+01, 6.7320000000e+01, 6.7390000000e+01, 6.7350000000e+01,
6.7360000000e+01, 6.7570000000e+01, 6.7660000000e+01, 6.7640000000e+01, 6.7800000000e+01, 6.7600000000e+01, 6.7900000000e+01, 6.7870000000e+01,
6.7940000000e+01, 6.7820000000e+01, 6.7950000000e+01, 6.8540000000e+01
... 500343 more values
} # data_jpeg2000_packing codedValues
====================== SECTION_8 ( length=4, padding=0 ) ======================
1-4 7777 = 7777从上述grib_dump的输出可以看到，Section 6中的bitMapIndicator被设置为0，表示使用位图段。
bitmap中保存了129600 bits的数据。
codedValues 中保存了500443个值，说明bitmap中应该由500443个数值为1，其余都是0。
下面使用简单的python代码来验证一下。
测试
使用perillaroc/nuwe-pyeccodes提供的 Python 3 接口读取 GRIB2 文件。
使用nuwe-pyeccodes获取message handler。
file_handler = nuwe_pyeccodes.GribFileHandler()
file_handler.openFile("bitmap.grib2")
message_handler = file_handler.next()读取bitmap数组，转型为np.bool，并打印值为true的个数，即数据点的个数。
bitmap_bytes = message_handler.getDoubleArray("bitmap")
bitmap_bytes = bitmap_bytes.astype(np.bool)
print(np.sum(bitmap_bytes == true)) # print 500443读取编码后的codedValues数组，并打印个数。
coded_values = message_handler.getDoubleArray("codedValues")
print(len(coded_values)) # print 500443将codesValues填入实际矩阵中，并打印图像。
values = np.zeros(1440*720)
values[bitmap_bytes] = coded_values
values = values.reshape(720, 1440)
imshow(values)
colorbar()
show()图片结果如下：
示例要素场绘图结果
使用numpy解析位图
上节中使用eccodes解码bitmap，下面使用python手动解码位图。
f = open("bitmap.grib2", "rb")
f.seek(16 + 21 + 72 + 34 + 23 + 6)
section_6 = f.read(129600)
unit_map_array = np.frombuffer(section_6, dtype=np.uint8)
bit_map_array = np.unpackbits(unit_map_array)bit_map_array与上节中的bitmap_bytes一致。
