GRIB2からマップタイルを出力したい

この記事は[防災アプリ開発 Advent Calendar 2024](https://adventar.org/calendars/9939)の20日目です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/311703/9d9f02a2-8935-15ef-ab1e-ed312ffc646c.png)
（2024年09月21日11時00分、全国合成レーダー降水強度）

## はじめに
普段目にする気象警報・注意報、地震情報や天気予報、雨雲レーダー、危険分布等のデータは、気象庁（気象業務支援センターが仲介）からがリアルタイムに配信されていますが、そのデータは直接人が見たり読めたりするものではありません。
データの容量が大きくなる雨雲レーダー、危険分布等のGPV（Grid Point Value)は特にそうです。

このようなGPVデータはマップに表示して見たいわけですが、調べてもデータをマップタイル（ラスタタイル）で出力してくれるプログラム等や知見、解説記事等は見当たりませんでした。
なので、一からやる必要があります。

## GRIB2とは

気象庁から配信されるGPV系のデータはほとんどGRIB2で配信されています。

GRIB2自体の説明はうまいことできる気がしないので、詳しい方の詳しい記事を紹介します。

https://qiita.com/e_toyoda/items/ce7497e1a633b16f1ff1

ちなみに、気象業務支援センターから配信されているデータは結構お高いです（データの契約数等によりけりだが、年間数万円から数百万円）。

## ラスタタイルを

色々試行錯誤する中で、以下の仕様で出力することにしました。

* TypeScript（Node.jsで実行）で構築
* 1タイルはタテ:512px、ヨコ:512px
* ズームレベルは、4から10（250mメッシュを表現するため、最大のズームレベルは最低でも10は必要）
* ズームレベル10で計算・作図し、9から4は一つ大きいズームレベルのタイルデータの色をもらって縮小
* 出力する画像はPNG
* PNGのインデックスカラーで容量削減、計算量削減


### じっそう

Node.jsの標準ライブラリだけで動作するよう、タイル、PNGの生成等はほかのパッケージを使わないように実装します。

高解像度降水ナウキャストや全国合成レーダー降水強度のように、格子のサイズが複数（陸域の250mと海域の1km）混在するデータの場合、同一現象、時間（実況・予想）でも第3節から第7節が複数回出力されるのでそれをうまく統一する処理、
また地理データを使いまわす（繰り返しの2回目以降の第3節を省略）して複数のデータを格納する場合もあることを考慮に入れる必要があります。

ラスタタイルを出力するPNGには、インデックスカラー（イメージデータ、8bitの場合は1byteで1px）で表現するので256色以上は使えないのであまり色が多くならないように間引く必要がありますが、今回はそんなデータないものとして扱います。

ビットマップの処理は省略します。（今回の場合では必要ない）


※最適化処理はとても長くなるので省略しています。

<details>
<summary>ソースコード</summary>

Node.js v22.2.0, v20.15.0 以降

```typescript:grib2.ts
import { EventEmitter } from 'node:events';
import { Readable } from 'node:stream';
import { MapTile } from './maptile';

export class GRIB2 extends EventEmitter<EventMap> {
  #buffer;
  #option;
  #dataLength = 0n;
  #header?: {
    center: number;
    subCenter: number;
    masterTableVersion: number;
    localTableVersion: number;
    dateTimeType: number;
    dateTime: Date;
    status: number;
    infoKind: number
  };
  #repetitionData?: {
    endTime: Date;
    startTime: Date;
    optional?: Record<string, any>;
    parameterCategory: number;
    parameterKind: number;
    processType: number;
    fixedSurface1: {
      kind: number | null;
      value: number | null;
    };
    fixedSurface2: {
      kind: number | null;
      value: number | null;
    };
    statisticalType: number | null;
    mapTile: MapTile;
    config: GRIB2OptionItem;
  };
  #section3?: Buffer;

  constructor(buffer: Buffer, option: GRIB2Option) {
    super();

    this.#buffer = buffer;
    this.#option = option;

    this.#sectionRead();
  }

  async #sectionRead() {
    console.time('MAIN');
    console.time('head read');
    this.#section0(await this.#buffer.subarray(0, 16));
    this.#section1(await this.#buffer.subarray(16, 37));
    console.timeEnd('head read');
    console.time('section read');

    for (const { section3, section4, section5, section6, section7 } of this.#repetition()) {
      await this.#sectionMain(section3, section4, section5, section6, section7);
    }

    console.timeEnd('section read');

    await this.#export();

    console.timeEnd('MAIN');
  }


  async #sectionMain(section3: Buffer, section4: Buffer, section5: Buffer, section6: Buffer, section7: Buffer) {
    const productTemplateKind = section4.readUintBE(7, 2);
    const parameterCategory = section4.readUInt8(9);
    const parameterKind = section4.readUInt8(10);
    const processType = section4.readUInt8(11);

    const fixedSurface1Kind = section4.readUInt8(22) === 0xff ? null : section4.readUInt8(22);
    const fixedSurface1ValueScale = section4.readUInt8(23);
    const fixedSurface1ValueFactor = section4.readUIntBE(24, 4);
    const fixedSurface1Value = (fixedSurface1ValueScale === 0xff || fixedSurface1ValueFactor === 0xffff_ffff) ? null :
      fixedSurface1ValueFactor / (10 ** fixedSurface1ValueScale);

    const fixedSurface2Kind = section4.readUInt8(28) === 0xff ? null : section4.readUInt8(28);
    const fixedSurface2ValueScale = section4.readUInt8(29);
    const fixedSurface2ValueFactor = section4.readUIntBE(30, 4);
    const fixedSurface2Value = (fixedSurface2ValueScale === 0xff || fixedSurface2ValueFactor === 0xffff_ffff) ? null :
      fixedSurface2ValueFactor / (10 ** fixedSurface2ValueScale);

    const statisticalType = productTemplateKind === 8 ? section4.readUInt8(46) : null;

    const targetConfig = this.#option.items.find(({ filter }) =>
      filter.category === parameterCategory && filter.kind === parameterKind && filter.processType === processType &&
      (!filter.fixedSurface1 || (filter.fixedSurface1.kind === fixedSurface1Kind && filter.fixedSurface1.value === fixedSurface1Value)) &&
      (!filter.fixedSurface2 || (filter.fixedSurface2.kind === fixedSurface2Kind && filter.fixedSurface2.value === fixedSurface2Value)) &&
      (typeof filter.statisticalType !== 'number' || filter.statisticalType === statisticalType)
    );

    if (!targetConfig || !this.#header) {
      return;
    }


    const {
      startTime,
      endTime,
      optional
    } = templateSection4(productTemplateKind, section4.subarray(9), this.#header.dateTime) ?? { endTime: new Date() };


    const data = this.#repetitionData;
    if (!data ||
      data.endTime.toJSON() !== endTime.toJSON() ||
      data.parameterCategory !== parameterCategory ||
      data.parameterKind !== parameterKind ||
      data.processType !== processType ||
      data.fixedSurface1.kind !== fixedSurface1Kind ||
      data.fixedSurface1.value !== fixedSurface1Value ||
      data.fixedSurface2.kind !== fixedSurface2Kind ||
      data.fixedSurface2.value !== fixedSurface2Value ||
      data.statisticalType !== statisticalType
    ) {
      await this.#export();
    }

    const mapTile = this.#repetitionData?.mapTile ??
      new MapTile({ ...this.#option, tileSize: 512 });

    this.#repetitionData ??= {
      startTime: startTime ?? endTime,
      endTime,
      optional,
      parameterCategory,
      parameterKind,
      processType,
      fixedSurface1: {
        kind: fixedSurface1Kind,
        value: fixedSurface1Value
      },
      fixedSurface2: {
        kind: fixedSurface2Kind,
        value: fixedSurface2Value
      },
      statisticalType,
      mapTile,
      config: targetConfig
    };

    const gridPointLonCount = section3.readUIntBE(30, 4);
    const gridPointLatCount = section3.readUIntBE(34, 4);
    const startLatitude = section3.readUIntBE(46, 4);
    const startLongitude = section3.readUIntBE(50, 4);
    const endLatitude = section3.readUIntBE(55, 4);
    const endLongitude = section3.readUIntBE(59, 4);
    const addLatitude = (startLatitude > endLatitude ? startLatitude - endLatitude : endLatitude - startLatitude) / (gridPointLatCount - 1) * (startLatitude > endLatitude ? 1 : -1);
    const addLongitude = (startLongitude > endLongitude ? startLongitude - endLongitude : endLongitude - startLongitude) / (gridPointLonCount - 1) * (startLongitude > endLongitude ? -1 : 1);

    const encoding = section5.readUIntBE(9, 2);

    if (encoding === 0) {
      // 単純圧縮
    }
    if (encoding === 200) {
      // ランレングス圧縮
      this.#RLE(
        startLongitude,
        startLatitude,
        addLongitude,
        addLatitude,
        gridPointLonCount,
        targetConfig,
        section5,
        section6,
        section7
      );
    }
  }

  #RLE(startLongitude: number, startLatitude: number, addLongitude: number, addLatitude: number, gridPointLonCount: number, config: GRIB2OptionItem, section5: Buffer, section6: Buffer, section7: Buffer) {
    const nBit = section5.readUInt8(11);
    const maxValue = section5.readUInt16BE(12);
    const M = section5.readUInt16BE(14);
    const lngu = 2 ** nBit - 1 - maxValue;
    const levelMap: { value: number; original: number; }[] = [{ value: 0, original: 0 }];
    const scalingFactor = 10 ** section5.readUInt8(16);

    for (let i = 17; i < section5.length; i += 2) {
      const v = section5.readUInt16BE(i);
      const original = (v & 0x7fff) * (v >> 15 === 1 ? -1 : 1);

      levelMap.push({ value: original / scalingFactor, original });
    }

    const colorMap: (string | null)[] = [null];
    const colorValues: [string, number][] = [];

    const colors = { type: 'level', min: 0, max: M };

    for (let i = 1; i < levelMap.length; i++) {
      const color = colorMap[i] = 
        (levelMap[i].value > 0 || this.#repetitionData?.config.isZeroValueColor) && levelMap[i].value >= colors.min ?
          colorValueRGB(((colors.type === 'value' ? levelMap[i].value : i) - colors.min) / (colors.max - colors.min)) : null;

      if (color) {
        colorValues.push([color, levelMap[i].original]);
      }
    }

    this.#repetitionData?.mapTile.setColors(colorValues, scalingFactor);

    const rle: [number, number, number][] = [];

    const horn1 = addLongitude / 2;
    const horn2 = addLatitude / 2;
    
    for (let i = 5, vn = 0, color = null as string | null, length = 0, count = 0, longitude = startLongitude, latitude = startLatitude; i < section7.length; i++) {
      const value = section7[i];
      if (value > maxValue) {
        length = lngu ** vn * (value - maxValue - 1);
        vn += 1;
      } else {
        color = colorMap[value];
        length = 1;
        vn = 0;
      }
      
      if (!color) {
        count += length;

        latitude = latitude - ~~(count / gridPointLonCount) * addLatitude;
        count = count % gridPointLonCount;

        longitude = startLongitude + count * addLongitude;
        continue;
      }

      while (length > 0) {
        const lineRemainingLength = (gridPointLonCount - count);
        const startRemainingLongitude = longitude;

        if (length >= lineRemainingLength) {
          longitude += addLongitude * (lineRemainingLength - 1);
          length -= lineRemainingLength;
          count += lineRemainingLength;
        } else {
          longitude += addLongitude * (length - 1);
          count += length;
          length = 0;
        }

        this.#repetitionData?.mapTile.rectFill(
          [(startRemainingLongitude - horn1) / 10e5, (latitude + horn2) / 10e5],
          [(longitude + horn1) / 10e5, (latitude - horn2) / 10e5],
          color
        );

        longitude += addLongitude;
        if (count >= gridPointLonCount) {
          latitude -= addLatitude;
          longitude = startLongitude;
          count = 0;
        }
      }
    }
  }



  async #export() {
    const header = this.#header;
    const data = this.#repetitionData;
    if (!header || !data) {
      return;
    }

    const tiles = await data.mapTile.export();
    
    const meta = {
      config: data.config,
      targetTime: header.dateTime,
      startTime: data.startTime,
      endTime: data.endTime,
      optional: data.optional
    };
    
    this.emit('export', meta , tiles);

    this.#repetitionData = void 0;
  }

  #section0(buffer: Buffer) {
    const grib = buffer.subarray(0, 4).toString('utf8');
    const gribVersion = buffer.readUInt8(7);

    if (grib !== 'GRIB' || gribVersion !== 2) {
      throw new Error('Not grib2 data stream.');
    }

    this.#dataLength = buffer.readBigUInt64BE(8);
  }

  #section1(buffer: Buffer) {
    const center = buffer.readUIntBE(5, 2);
    const subCenter = buffer.readUIntBE(7, 2);
    const masterTableVersion = buffer.readUInt8(9);
    const localTableVersion = buffer.readUInt8(10);
    const dateTimeType = buffer.readUInt8(11);

    const dataTime = new Date(
      buffer.readUint16BE(12),
      buffer.readUint8(14) - 1,
      buffer.readUint8(15),
      buffer.readUint8(16),
      buffer.readUint8(17),
      buffer.readUint8(18),
      0
    );
    dataTime.setMinutes(dataTime.getMinutes() - dataTime.getTimezoneOffset());

    const status = buffer.readUInt8(19);
    const infoKind = buffer.readUInt8(20);

    this.#header = {
      center,
      subCenter,
      masterTableVersion,
      localTableVersion,
      dateTime: dataTime,
      dateTimeType,
      status,
      infoKind
    };
  }

  #repetition() {
    const list:{ section3: Buffer; section4: Buffer; section5: Buffer; section6: Buffer; section7: Buffer; }[] = [];

    let start = 37;
    while (true) {
      const section3head = this.#buffer.subarray(start, start + 5);
      const section3Length = section3head.readUInt32BE();

      // END SECTION
      if (section3Length === 926365495) {
        break;
      }


      const section3 = section3head.readUInt8(4) === 3 ? this.#buffer.subarray(start, start += section3Length) : this.#section3;

      if (!section3 || section3.readUInt8(4) === 2) {
        continue;
      }

      this.#section3 = section3;

      const section4Length = this.#buffer.subarray(start, start + 4).readUInt32BE();

      const section4 = this.#buffer.subarray(start, start += section4Length);

      const section5Length = this.#buffer.subarray(start, start + 4).readUInt32BE();

      const section5 = this.#buffer.subarray(start, start += section5Length);

      const section6Length = this.#buffer.subarray(start, start + 4).readUInt32BE();

      const section6 = this.#buffer.subarray(start, start += section6Length);

      const section7Length = this.#buffer.subarray(start, start + 4).readUInt32BE();

      const section7 = this.#buffer.subarray(start, start += section7Length);

      list.push({ section3, section4, section5, section6, section7 });
    }

    return list;
  }
}

function colorValueRGB(value: number) {
  const rgb = [0, 0, 0];

  if (value >= 1) {
    rgb[0] = 255;
    rgb[1] = 0;
    rgb[2] = 0;
  } else if (value >= 0.75) {
    rgb[0] = 255;
    rgb[1] = ~~((1 - (value - 0.75) / 0.25) * 255);
    rgb[2] = 0;
  } else if (value >= 0.5) {
    rgb[0] = ~~((value - 0.5) / 0.25 * 255);
    rgb[1] = 255;
    rgb[2] = 0;
  } else if (value >= 0.25) {
    rgb[0] = 0;
    rgb[1] = 255;
    rgb[2] = ~~((1 - (value - 0.25) / 0.25) * 255);
  } else if (value >= 0) {
    rgb[0] = 0;
    rgb[1] = ~~(value / 0.25 * 255);
    rgb[2] = 255;
  } else {
    rgb[0] = 0;
    rgb[1] = 0;
    rgb[2] = 255;
  }

  return `#${rgb.map(v => v.toString(16).padStart(2, '0')).join('')}`;
}

function templateSection4(num: number, data: Buffer, targetTime = new Date()): { startTime?: 
 Date; endTime: Date; optional?: any } | null {
  if (num === 0) {
    return templateSection4_0(data, targetTime);
  }
  if (num === 8 || num === 50000 || num === 50008 || num === 50009 || num === 50010 || num === 50011 || num === 50012) {
    return templateSection4_8(data, targetTime);
  }
  return null;
}

function templateSection4_0(data: Buffer, targetTime: Date) {
  const forecastType = data.readUInt8(8);
  const forecastValue = data.readUIntBE(9, 4);

  const endTime = new Date(targetTime);

  const time = (forecastValue >>> 31) === 1 ? -(forecastValue & 0x7fffffff) : value;

  if (forecastType === 0) {
    resultDate.setMinutes(resultDate.getMinutes() + time);
  }
  if (forecastType === 1) {
    resultDate.setHours(resultDate.getHours() + time);
  }
  if (forecastType === 2) {
    resultDate.setDate(resultDate.getDate() + time);
  }
  if (forecastType === 3) {
    resultDate.setMonth(resultDate.getMonth() + time);
  }
  if (forecastType === 4) {
    resultDate.setFullYear(resultDate.getFullYear() + time);
  }
  if (forecastType === 5) {
    resultDate.setFullYear(resultDate.getFullYear() + time * 10);
  }
  if (forecastType === 6) {
    resultDate.setFullYear(resultDate.getFullYear() + time * 30);
  }
  if (forecastType === 7) {
    resultDate.setFullYear(resultDate.getFullYear() + time * 100);
  }
  if (forecastType === 10) {
    resultDate.setHours(resultDate.getHours() + time * 3);
  }
  if (forecastType === 11) {
    resultDate.setHours(resultDate.getHours() + time * 6);
  }
  if (forecastType === 12) {
    resultDate.setHours(resultDate.getHours() + time * 12);
  }
  if (forecastType === 13) {
    resultDate.setSeconds(resultDate.getSeconds() + time);
  }

  return { endTime };
}

function templateSection4_8(data: Buffer, targetTime: Date) {
  const { endTime: startTime } = templateSection4_0(data, targetTime);

  const endTime = new Date();

  endTime.setFullYear(
    data.readUint16BE(25),
    data.readUint8(27) - 1,
    data.readUint8(28)
  );
  endTime.setHours(
    data.readUint8(29),
    data.readUint8(30),
    data.readUint8(31),
    0
  );

  endTime.setMinutes(endTime.getMinutes() - endTime.getTimezoneOffset());

  return { startTime, endTime };
}

export interface GRIB2Option {
  maxZoom: number;
  minZoom: number;
  bounds: [[number, number], [number, number]];
  items: GRIB2OptionItem[];
}

interface GRIB2OptionItem {
  filter: {
    category: number;
    kind: number;
    processType: number;
    fixedSurface1?: {
      kind: number;
      value: number | null;
    };
    fixedSurface2?: {
      kind: number;
      value: number | null;
    };
    statisticalType?: number;
  };
  isZeroValueColor?: boolean;
}

interface EventMap {
  export: [
    {
      config: GRIB2OptionItem;
      targetTime: Date;
      startTime: Date;
      endTime: Date;
      optional?: Record<string, any>;
    },
    [`${number}_${number}_${number}` ,Buffer][]
  ];
}
```

```typescript:maptile.ts
import { PNG } from './png';

export class MapTile {
  #option;
  #tileSize;
  #bounds;
  #zoom;
  #pngs = new Map<`${number}_${number}`, PNG>();
  #pngColorValues: [string, number][] = [];
  #pngColors: string[] = [];
  #pngColorValueScalingFactor = 1;

  constructor(option: MapTileOption) {
    if (typeof option.zoom === 'number' && (option.zoom > option.maxZoom || option.zoom < option.minZoom)) {
      throw new Error(`option error: [zoom: ${option.zoom}] range: ${option.minZoom} ~ ${option.maxZoom}`);
    }

    this.#option = option;
    this.#tileSize = option.tileSize ?? 256;
    this.#zoom = option.zoom ?? option.maxZoom;

    this.#bounds = {
      start: [
        Math.min(...option.bounds.map(([lon, lat]) => lon)),
        Math.min(...option.bounds.map(([lon, lat]) => lat))
      ],
      end: [
        Math.max(...option.bounds.map(([lon, lat]) => lon)),
        Math.max(...option.bounds.map(([lon, lat]) => lat))
      ]
    };
  }

  setColors(colors: [string, number][], scalingFactor = 0) {
    this.#pngColorValues = colors;
    this.#pngColors = colors.map(([color]) => color);
    this.#pngColorValueScalingFactor = scalingFactor;
  }

  rectFill(startLocation: [number, number], endLocation: [number, number], color: string) {
    const [startLon, startLat] = startLocation;
    const [endLon, endLat] = endLocation;

    if (startLon < this.#bounds.start[0] || startLon > this.#bounds.end[0] ||
      startLat < this.#bounds.start[1] || startLat > this.#bounds.end[1] ||
      endLon < this.#bounds.start[0] || endLon > this.#bounds.end[0] ||
      endLat < this.#bounds.start[1] || endLat > this.#bounds.end[1]) {
      return;
    }

    const startPointX = lon2point(
      startLon < endLon ? startLon : endLon,
      this.#zoom
    );
    const startPointY = lat2point(
      startLat < endLat ? endLat : startLat,
      this.#zoom
    );
    const endPointX = lon2point(
      startLon < endLon ? endLon : startLon,
      this.#zoom
    );
    const endPointY = lat2point(
      startLat < endLat ? startLat : endLat,
      this.#zoom
    );

    const startTileNumX = ~~(startPointX / this.#tileSize);
    const startTileNumY = ~~(startPointY / this.#tileSize);
    const endTileNumX = endPointX % this.#tileSize === 0 ? ~~(endPointX / this.#tileSize) - 1 : ~~(endPointX / this.#tileSize);
    const endTileNumY = endPointY % this.#tileSize === 0 ? ~~(endPointY / this.#tileSize) - 1 : ~~(endPointY / this.#tileSize);

    const x = startPointX % this.#tileSize;
    const y = startPointY % this.#tileSize;
    const w = endPointX - startPointX;
    const h = endPointY - startPointY;


    for (let tileX = startTileNumX; tileX <= endTileNumX; tileX++) {
      for (let tileY = startTileNumY; tileY <= endTileNumY; tileY++) {
        const dX = (startTileNumX - tileX) * this.#tileSize;
        const dY = (startTileNumY - tileY) * this.#tileSize;

        const key = `${tileX}_${tileY}` as const;

        const png = this.#pngs.get(key) ??
          this.#pngs.set(key, new PNG(this.#tileSize, this.#tileSize, this.#pngColors)).get(key);

        png?.rectFill(
          dX + x,
          dY + y,
          w,
          h,
          color
        );
      }
    }
  }

  async export() {
    return await this.#export();
  }

  async #export(): Promise<[`${number}_${number}_${number}`, Buffer][]> {
    const outputBuffersPromise = this.#pngBuffers();
    if (this.#zoom <= this.#option.minZoom) {
      return await outputBuffersPromise;
    }

    const buffers = await Promise.all([
      outputBuffersPromise,
      this.#pngHalf().then(() => this.#export())
    ]);


    return buffers.flat();
  }

  async #pngBuffers() {
    const zoom = this.#zoom;

    const pngs = [...this.#pngs];

    const outputBuffers: [`${number}_${number}_${number}`, Buffer][] = [];

    for (let i = 0; i < pngs.length; i++) {
      const [tileName, png] = pngs[i];
      const buffer = await png.export();

      if (!buffer){
        continue;
      }

      outputBuffers.push([
        `${zoom}_${tileName}`,
        buffer   
      ]);
    }

    return outputBuffers;
  }

  async #pngHalf() {
    console.time('halfsize');

    const pngs = [...this.#pngs];

    const pngMap = new Map<`${number}_${number}`, PNG>();

    for (let i = 0; i < pngs.length; i++) {
      const [tileName, png] = pngs[i];

      const [oX, oY] = tileName.split('_').map(v => +v);

      const cX = ~~(oX / 2);
      const cY = ~~(oY / 2);

      const key = `${cX}_${cY}` as const;

      const hPng = pngMap.get(key) ??
        pngMap.set(key, new PNG(this.#tileSize, this.#tileSize)).get(key);

      const size = this.#tileSize / 2;

      hPng?.copyHalfSize(png, (oX % 2) * size, (oY % 2) * size);
    }

    pngs.length = 0;
    this.#zoom = this.#zoom - 1;

    this.#pngs = pngMap;

    console.timeEnd('halfsize');
  }
}


export interface MapTileOption {
  tileSize?: number;
  bounds: [[number, number], [number, number]];
  minZoom: number;
  maxZoom: number;
  zoom?: number;
}

function lat2point(lat: number, zoom: number) {
  const p = Math.sin(Math.PI / 180 * lat);
  const z = 2 ** (zoom + 7) / Math.PI;

  return Math.round(z * (-(Math.log((1 + p) / (1 - p)) / 2) + Math.atanh(Math.sin(Math.PI / 180 * 85.05112878))));
}

function lon2point(lon: number, zoom: number) {
  return Math.round(2 ** (zoom + 7) * (lon / 180 + 1));
}

```
```typescript:png.ts
import { promisify } from 'node:util';
import { deflate, crc32 } from 'node:zlib';

const deflatePromise = promisify(deflate);

export class PNG {
  #width;
  #height;
  #buffer;
  #colors;
  #colorsIndex: { [key: string]: number } = {};

  constructor(width: number, height: number, colors: string[] = [], background = '#ffffff00') {
    this.#width = ~~width;
    this.#height = ~~height;

    this.#buffer = Buffer.alloc((this.#width + 1) * this.#height);
    this.#colors = [background, ...colors];

    for (let i = 0; i < this.#colors.length; i++) {
      this.#colorsIndex[this.#colors[i]] = i;
    }
  }

  rectFill(x: number, y: number, width: number, height: number, color: string) {
    const fillWidth = ~~(
      x >= 0 ?
        (x + width > this.#width ? this.#width - x : width) :
        (x + width > this.#width ? this.#width : x + width)
    );
    const fillHeight = ~~(
      y >= 0 ?
        (y + height > this.#height ? this.#height - y : height) :
        (y + height > this.#height ? this.#height : y + height)
    );

    if (x < 0) {
      x = 0;
    }
    if (y < 0) {
      y = 0;
    }
    x = ~~x;
    y = ~~y;

    let colorIndex = this.#colorsIndex[color];
    if (!colorIndex) {
      this.#colors.push(color);
      colorIndex = this.#colors.indexOf(color);
      this.#colorsIndex[color] = colorIndex;

      if (colorIndex > 255) {
        throw new Error('color index overflow');
      }
    }

    for (let i = 0; i < fillHeight; i++) {
      const offset = (i + y) * (this.#width + 1) + x + 1;
      const count = offset + fillWidth;

      for (let j = offset; j < count; j++) {
        this.#buffer[j] = colorIndex;
      }
    }
  }

  copyHalfSize(src: PNG, x: number, y: number) {
    const srcWidthHalf = src.#width / 2;
    const srcHeightHalf = src.#height / 2;

    if (src.#width !== src.#height || this.#width < x + srcWidthHalf || this.#height < y + srcHeightHalf) {
      return;
    }

    for (let i = 0; i < src.#colors.length; i++) {
      const color = src.#colors[i];

      if (typeof this.#colorsIndex[color] !== 'number') {
        this.#colors.push(color);
        this.#colorsIndex[color] = this.#colors.indexOf(color);
      }
    }

    for (let sHY = 0; sHY < srcHeightHalf; sHY++) {
      const sY = sHY * 2 * (src.#width + 1);
      for (let sHX = 0; sHX < srcWidthHalf; sHX++) {
        const sX = sHX * 2 + 1;

        const value = src.#buffer[sY + sX];

        if (value === 0) {
          continue;
        }

        const color = src.#colors[value];

        const colorIndex = this.#colorsIndex[color];

        const bi = x + sHX + (y + sHY) * (this.#width + 1) + 1;
        this.#buffer[bi] = colorIndex;
      }
    }

  }

  async export() {
    if (this.#colors.length <= 1) {
      return;
    }

    const IHDR = Buffer.alloc(25);
    IHDR.writeUIntBE(13, 0, 4);
    IHDR.write('IHDR', 4, 4);
    IHDR.writeUIntBE(this.#width, 8, 4);
    IHDR.writeUIntBE(this.#height, 12, 4);
    IHDR.writeUInt8(8, 16);
    IHDR.writeUInt8(3, 17);
    IHDR.writeUInt8(0, 18);
    IHDR.writeUInt8(0, 19);
    IHDR.writeUInt8(0, 20);

    const PLTE = Buffer.alloc((this.#colors.length * 3) + 12);
    PLTE.writeUIntBE(this.#colors.length * 3, 0, 4);
    PLTE.write('PLTE', 4, 4);

    const tRNS = Buffer.alloc((this.#colors.length) + 12);
    tRNS.writeUIntBE(this.#colors.length, 0, 4);
    tRNS.write('tRNS', 4, 4);

    for (let i = 0; i < this.#colors.length; i++) {
      const [, r, g, b, a] = this.#colors[i]
        .match(/^#([a-z0-9]{2})([a-z0-9]{2})([a-z0-9]{2})([a-z0-9]{2})?$/i) ?? [];

      PLTE.set(
        [
          parseInt(r ?? 0, 16),
          parseInt(g ?? 0, 16),
          parseInt(b ?? 0, 16)
        ],
        i * 3 + 8);

      tRNS.writeUInt8(parseInt(a ?? 'ff', 16), i + 8);
    }

    const data = await deflatePromise(this.#buffer, { level: 5 });

    const IDAT = Buffer.alloc(data.length + 12);
    IDAT.writeUIntBE(data.length, 0, 4);
    IDAT.write('IDAT', 4, 4);
    IDAT.set(data, 8);

    const buffers = [
      IHDR,
      PLTE,
      tRNS,
      IDAT
    ];

    for (const buffer of buffers) {
      buffer.writeUIntBE(crc32(buffer.subarray(4, buffer.length - 4)), buffer.length - 4, 4);
    }

    return Buffer.concat([
      Buffer.from([0x89, 0x50, 0x4e, 0x47, 0x0d, 0xa, 0x1a, 0xa]),
      ...buffers,
      Buffer.from([0, 0, 0, 0, 0x49, 0x45, 0x4e, 0x44, 0xae, 0x42, 0x60, 0x82])
    ]);
  }
}
```
</details>

----
### 実行と出力されたデータ

```typescript
import { readFileSync } from 'node:fs';
import { GRIB2 } from './grib2';

// 高解像度降水ナウキャスト 降水強度 実況値を抜き出し
new GRIB2(readFileSync('./Z__C_RJTD_yyyyMMddhhmmss_NOWC_GPV_Ggis0p25km_Pri60lv_Aper5min_FH0000-0030_grib2.bin'), {
  bounds: [[118, 20], [150, 48]], // 描写範囲
  maxZoom: 10, // 出力するタイルの最大ズーム
  minZoom: 4, // 出力するタイルの最小ズーム
  items: [{ // 出力するデータの定義等
    filter: {  // データの絞り込み
      category: 1, // 第4節 パラメーターカテゴリー
      kind: 203, // 第4節 パラメーター番号
      processType: 0 // 第4節 作成処理の種類
    }
  }]
})
  .on('export', (meta, tiles) => console.log(meta, tiles))
```

気象庁のHPに[GPVサンプルデータの一覧](https://www.data.jma.go.jp/developer/gpv_sample.html)という便利なものがあるのでこれを使って試してみてください（GSM、MSM、LFMなどは単純圧縮データなので使えません）。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/311703/d1d4ec80-2da6-dfa5-e505-e92367caa8c1.png)
全国合成レーダー画像を出力したものをそのまま地図に乗せるとこうなり


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/311703/bde5844d-db57-2942-d860-d6b45bb26378.png)
ブラウザで気象庁カラーに変換するとこうなる

（本記事の地図データは OpenStreetMap のものを用いた）


----
### 苦労のタイムライン

(半年も前なので記憶はあやふや)

#### 2024/06/21

https://x.com/soshi1822/status/1804004842282848496

GRIB2のデータには無効域、降水量0などのデータも入っているわけですが、そのいらないデータを読み飛ばす処理を間違っていた時の図です。

https://x.com/soshi1822/status/1804018973132099890

ラスタタイル自体を生成するとき、タイル1枚につき1つのPNGを作成しそこにデータを書き込んでいくわけですが、タイルを複数に跨ぐデータが来た場合の処理が少し面倒だった記憶があります。

https://x.com/soshi1822/status/1804075587490975957

全国合成レーダー降水強度の場合、250mと1kmメッシュが混在する都合上いくつかの領域に分かれて第3節～第7節が入っています。
その領域境界のところに隙間が入ってしまいましたが、これは第3節の格子間隔を利用したためです。
格子系のメッシュサイズは緯度経度で入っているのですが割り切れない数が入っていることがあるため素直に使うと誤差が蓄積して隙間が発生しました。
最初の格子点と最後の格子点の緯度経度から格子のサイズを求め解決しています。

#### 2024/06/22

https://x.com/soshi1822/status/1804482196537045497

単純圧縮で配信されている（本記事では処理は書いていない）、北西太平洋高解像度日別海面水温解析格子点資料を解析していて、変なデータで困っているところです。
原因は単純なミスで、読み取り位置が本来の5byteほど前だったことです。

本来はこんな感じ↓

https://x.com/soshi1822/status/1804720094867841229

#### 2024/06/23

https://x.com/soshi1822/status/1804742889295262159

調子に乗ってGRIB2で配信されていない推計震度分布を出力したところです。
推計震度分布には最適化されていないので、格子点の数はそこまでではないものの処理時間がかかっているのを改善するのが今後の課題です。

#### 2024/08 ~

大体の実装を終え、鮮やかなグラデーションで出力された細かい値を持った元データをブラウザで荒いデータに変換する実装を行うのですが、長くなってしまうのでまたの機会にさせていただきたく思います。

まぁ、タイル作成時点で荒くさせてしまうのが一番手っ取り早いですが、いろんなスケールで見たい、基準値の超過率、平年との比較、格子の代表値を確認したいのでそのままの細かいままで作成しています。

以下はそんな例

https://x.com/soshi1822/status/1823743421154255039

天気予報分布図の最高気温で25度以上、30度以上は1度ごとに塗り分けたもの

https://x.com/soshi1822/status/1828806297283629499

大雨特別警報の基準超過率（土壌雨量指数）
基準がない、5％未満のところは透明です

## おわりに

せっかく[いろいろなデータを買って受信している](https://dmdata.jp/docs/jmafile/)のだから何か（サービス等）はしたいのだがあまりモチベーションと意義を見出せず、足踏みをしております。
何かいい活用方法等のお考えがありましたら是非ご意見ください。

以上、ありがとうございました。