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一致。