GRIB学习笔记：GRIB2解码 - 处理JPEG2000压缩
目录
在之前的一些列文章中，我使用中国气象局数值预报中心GRAPES全球模式生成的GRIB2产品介绍GRIB2格式。接下来的一组文章会介绍如何解析GRAPES模式生成的GRIB2数据。
尽管已有wgrib2、ecCodes等业务领域广泛应用的GRIB2文件编解码工具，自己实际编写一个GRIB2解码器有助于我更深入地了解GRIB2格式。
简要研究过ecCodes及其前身grib api的源码，为了能解析任意GRIB2文件，大量代码用于处理GRIB2的表格和模板，整个项目的结构比较复杂。数值预报中心生成的GRIB2数据往往使用相同的表格和模板。对于手动编写解码程序，仅需要考虑一种情况。因此本系列文章仅关注固定表格、模板组合下的GRIB2数据解码。
本文主要介绍GRIB2消息解码唯一需要外部库的部分：JPEG2000解码。
介绍
数值预报中心的GRIB2数据都使用JPEG2000压缩，详情参见《GRIB学习笔记：样例文件分析 - Section 5》。本文使用上面文章的样例数据进行说明。
ecCodes中支持使用Jasper和OpenJPEG解码JPEG2000数据，本文选择OpenJPEG解码数据，并参考ecCodes中的代码。
数据获取
数据的描述放在Section 5，包括数据点个数，参考值，二进制因子和十进制因子等。
int binary_scale_factor = 0x0;
int decimal_scale_factor = 0x2;
float convertToFloat(uint32_t v) {
return *(reinterpret_cast<float*>(&v));
}
float reference_value = convertToFloat(0x46B1298E);我使用64位Windows 10的Visual Studio 2019编译程序，float是4字节浮点数，所以用float表示参考值。
压缩后的数据放在Section 7，根据压缩的数据起始位置和终止位置，读取数据。
auto start_pos = 0xB1;
auto raw_data_length = 0x98998 - 0xB1 + 1;
auto buf = new unsigned char[raw_data_length];
std::FILE* f = std::fopen(grib_file_path.c_str(), "rb");
std::fseek(f, start_pos, SEEK_SET);
std::fread(buf, 1, raw_data_length, f);
std::fclose(f);解码
使用ecCodes的代码，使用OpenJPEG解码。需要三个输入数据：压缩后的数据及字节长度，原始数据个数。
std::vector<double> decodeValues(unsigned char* buf, size_t raw_data_length, size_t data_count) {
int err = 0;
unsigned long mask;
std::vector<double> val;
opj_dparameters_t parameters = { 0, }; /* decompression parameters */
opj_stream_t* stream = NULL;
opj_memory_stream mstream;
opj_image_t* image = NULL;
opj_codec_t* codec = NULL;
opj_image_comp_t comp = { 0, };
/* set decoding parameters to default values */
opj_set_default_decoder_parameters(¶meters);
parameters.decod_format = 1; /* JP2_FMT */
/* get a decoder handle */
codec = opj_create_decompress(OPJ_CODEC_J2K);
/* catch events using our callbacks and give a local context */
opj_set_info_handler(codec, openjpeg_info, nullptr);
opj_set_warning_handler(codec, openjpeg_warning, nullptr);
opj_set_error_handler(codec, openjpeg_error, nullptr);
/* initialize our memory stream */
mstream.pData = buf;
mstream.dataSize = raw_data_length;
mstream.offset = 0;
/* open a byte stream from memory stream */
stream = opj_stream_create_default_memory_stream(&mstream, OPJ_STREAM_READ);
/* setup the decoder decoding parameters using user parameters */
if (!opj_setup_decoder(codec, ¶meters)) {
err = 2;
goto cleanup;
}
if (!opj_read_header(stream, codec, &image)) {
err = 3;
goto cleanup;
}
if (!opj_decode(codec, stream, image)) {
err = 4;
goto cleanup;
}
if (!(data_count <= image->comps[0].w * image->comps[0].h)) {
err = 5;
goto cleanup;
}
if ((image->numcomps != 1) || (image->x1 * image->y1) == 0) {
err = 6;
goto cleanup;
}
_ASSERT(image->comps[0].sgnd == 0);
_ASSERT(comp.prec <= sizeof(image->comps[0].data[0]) * 8 - 1); /* BR: -1 because I don't know what happens if the sign bit is set */
_ASSERT(image->comps[0].prec < sizeof(mask) * 8 - 1);
auto data = image->comps[0].data;
mask = (1 << image->comps[0].prec) - 1;
auto count = image->comps[0].w * image->comps[0].h;
val.resize(count);
for (auto i = 0; i < count; i++) {
auto v = data[i];
val[i] = v & mask;
}
if (!opj_end_decompress(codec, stream)) {
err = 7;
}
cleanup:
/* close the byte stream */
if (codec) opj_destroy_codec(codec);
if (stream) opj_stream_destroy(stream);
if (image) opj_image_destroy(image);
return val;
}还原数值
解码得到的数据只是打包值，还需要根据公式还原原始数据。
for (auto i = 0; i < data_count; i++) {
val[i] = (reference_value + val[i] * std::pow(2, binary_scale_factor)) / std::pow(10, decimal_scale_factor);
}示例
使用GRAPES GFS模式输出的样例文件测试，解码结果如下图所示：
OpenJPEG封装
ecCodes使用 https://groups.google.com/forum/#!topic/openjpeg/8cebr0u7JgY 中封装OpenJPEG的函数，用于解码内存数组中的JPEG压缩数据，而不是来自图片文件的数据。
#include <openjpeg.h>
typedef struct j2k_encode_helper {
size_t buffer_size;
long width;
long height;
long bits_per_value;
float compression;
long no_values;
const double* values;
double reference_value;
double divisor;
double decimal;
long jpeg_length;
unsigned char* jpeg_buffer;
} j2k_encode_helper;
/* OpenJPEG 2.1 version of grib_openjpeg_encoding.c */
/* opj_* Helper code from https://groups.google.com/forum/#!topic/openjpeg/8cebr0u7JgY */
/* These routines are added to use memory instead of a file for input and output */
/* struct need to treat memory as a stream */
typedef struct
{
OPJ_UINT8* pData; /* our data */
OPJ_SIZE_T dataSize; /* how big is our data */
OPJ_SIZE_T offset; /* where we are currently in our data */
j2k_encode_helper* helper;
} opj_memory_stream;
/* This will read from our memory to the buffer */
OPJ_SIZE_T opj_memory_stream_read(void* buffer, OPJ_SIZE_T nb_bytes, void* p_user_data)
{
opj_memory_stream* mstream = (opj_memory_stream*)p_user_data; /* Our data */
OPJ_SIZE_T nb_bytes_read = nb_bytes; /* Amount to move to buffer */
/* Check if the current offset is outside our data buffer */
if (mstream->offset >= mstream->dataSize)
return (OPJ_SIZE_T)-1;
/* Check if we are reading more than we have */
if (nb_bytes > (mstream->dataSize - mstream->offset))
nb_bytes_read = mstream->dataSize - mstream->offset;
memcpy(buffer, &(mstream->pData[mstream->offset]), nb_bytes_read);
mstream->offset += nb_bytes_read; /* Update the pointer to the new location */
return nb_bytes_read;
}
/* Write from the buffer to our memory */
OPJ_SIZE_T opj_memory_stream_write(void* buffer, OPJ_SIZE_T nb_bytes, void* user_data)
{
opj_memory_stream* mstream = (opj_memory_stream*)user_data; /* our data */
OPJ_SIZE_T nb_bytes_write = nb_bytes; /* Amount to move to buffer */
/* Check if the current offset is outside our data buffer */
if (mstream->offset >= mstream->dataSize)
return (OPJ_SIZE_T)-1;
/* Check if we are writing more than we have space for */
if (nb_bytes > (mstream->dataSize - mstream->offset))
nb_bytes_write = mstream->dataSize - mstream->offset;
/* Copy the data from the internal buffer */
memcpy(&(mstream->pData[mstream->offset]), buffer, nb_bytes_write);
mstream->offset += nb_bytes_write; /* Update the pointer to the new location */
return nb_bytes_write;
}
/* Moves the pointer forward, but never more than we have */
OPJ_OFF_T opj_memory_stream_skip(OPJ_OFF_T nb_bytes, void* user_data)
{
opj_memory_stream* mstream = (opj_memory_stream*)user_data;
OPJ_SIZE_T l_nb_bytes;
if (nb_bytes < 0)
return -1; /* No skipping backwards */
l_nb_bytes = (OPJ_SIZE_T)nb_bytes; /* Allowed because it is positive */
/* Do not allow jumping past the end */
if (l_nb_bytes > mstream->dataSize - mstream->offset)
nb_bytes = mstream->dataSize - mstream->offset;
mstream->offset += l_nb_bytes;
return l_nb_bytes; /* Returm how far we jumped */
}
/* Sets the pointer to anywhere in the memory */
OPJ_BOOL opj_memory_stream_seek(OPJ_OFF_T nb_bytes, void* user_data)
{
opj_memory_stream* mstream = (opj_memory_stream*)user_data;
if (nb_bytes < 0)
return OPJ_FALSE; /* Not before the buffer */
if (nb_bytes > (OPJ_OFF_T) mstream->dataSize)
return OPJ_FALSE; /* Not after the buffer */
mstream->offset = (OPJ_SIZE_T)nb_bytes; /* Move to new position */
return OPJ_TRUE;
}
void opj_memory_stream_do_nothing(void* p_user_data)
{
OPJ_ARG_NOT_USED(p_user_data);
}
/* Create a stream to use memory as the input or output */
opj_stream_t* opj_stream_create_default_memory_stream(opj_memory_stream* memoryStream, OPJ_BOOL is_read_stream)
{
opj_stream_t* stream;
if (!(stream = opj_stream_default_create(is_read_stream)))
return (NULL);
/* Set how to work with the frame buffer */
if (is_read_stream)
opj_stream_set_read_function(stream, opj_memory_stream_read);
else
opj_stream_set_write_function(stream, opj_memory_stream_write);
opj_stream_set_seek_function(stream, opj_memory_stream_seek);
opj_stream_set_skip_function(stream, opj_memory_stream_skip);
opj_stream_set_user_data(stream, memoryStream, opj_memory_stream_do_nothing);
opj_stream_set_user_data_length(stream, memoryStream->dataSize);
return stream;
}