Back to home

高解像度降水ナウキャストのGPVデータを読んでみる

8 min read
Table of Contents

お久しぶりです。2025年のお盆はいかがお過ごしでしたでしょうか。

私は、お客様から急遽依頼された短期間でのアプリ開発に取り組んでいます。

このような超特急での開発には、追加の報酬をいただきたいものです。

さて、愚痴はこれくらいにして、今回は気象庁が提供する高解像度降水ナウキャストのGPVデータを読み解いてみた話をご紹介します。

たまたま気象庁のGPVデータを扱う機会があったので、どのようなものか調べてみました。
※ 冒頭で触れたアプリ開発とは関係ありません。

注意

この記事は、私が調査した内容をまとめたものです。

Zennでいうスクラップ記事のようなものであり、一部に誤った情報が含まれている可能性があります。

もし間違いを見つけられた際は、X(旧Twitter)でお知らせいただけると幸いです。

高解像度降水ナウキャストとは

高解像度降水ナウキャストとは、降水強度および5分間降水量について、分布図形式で予報を行うものです。

250m格子単位で5分から30分先まで、1km格子単位で35分から60分先までの予報を提供しています。

降水ナウキャストというものもありますが、こちらも1時間後までの降水強度を分布図形式で予報するものです。

両者の違いは、提供している格子の種類です。

降水ナウキャストは1km格子単位のみですが、高解像度降水ナウキャストには250m格子と1km格子の両方が含まれているため、データを確認する際に非常に手間がかかります。

GPVデータとは

記事の冒頭で触れた GPV とは Grid Point Value の略称で、格子点上の気象要素などの値を表す気象データのことです。

詳細については、以下のリンクをご参照ください。

GPVから最も近い格子点を抽出する方法|お天気データサイエンス|日本気象株式会社
ods.n-kishou.co.jp
image
GPV形式の気象データ入門|加藤芳樹 (Yoshiki Kato)
note.com
image

ごく簡単に言えば、格子点ごとに気象データが格納されているものです。


出典: 気象庁

早速中身を見てみる

事前知識の説明はここまでにして、早速データの中身を見ていきましょう。

中身を見るための準備

GPVデータはGRIB2形式で提供されています。

そのため、GRIB2を読み込むためのツールが必要です。

ホストPCにGRIB2関連ツールを直接インストールすると余計なファイルが増えそうなので、Dockerコンテナを使用します。

Dockerコンテナの準備

Gitからクローンしてくる方法を採用します。

Dockerfileは以下のような構成です。
wgrib2 以外にも調査に必要なツールを含めています。

FROM python:3.12.3-slim

RUN apt-get update && apt-get install -y \
    build-essential \
    cmake \
    git \
    gzip \
    zip \
    unzip \
    libnetcdf-dev \
    libnetcdff-dev \
    netcdf-bin \
    gdal-bin \
    libgdal-dev \
    zlib1g-dev \
    libtiff-dev \
    && apt-get clean && rm -rf /var/lib/apt/lists/*

# Pythonパッケージのインストール
RUN pip install --no-cache-dir setuptools numpy netCDF4 rio-pmtiles


# Install wgrib2
RUN cd ~ && \
    git clone --depth 1 --branch v3.7.0 https://github.com/NOAA-EMC/wgrib2.git
RUN cd ~/wgrib2 \
    && mkdir build && cd build \
    && cmake .. -D USE_NETCDF=1 -D ENABLE_NETCDF=1 -D ENABLE_GRIB2=1 -D ENABLE_JASPER=0 \
    && make \
    && make install \
    && cp ./install/bin/wgrib2 /usr/local/bin \
    && cd ~ && rm -rf wgrib2

CMD ["tail", "-f", "/dev/null"]

docker-compose.yml は以下のようになります。

Docker Composeで起動する方が楽なので、こちらを使用します。

services:
  now_cast_amedas:
    build:
      context: ./docker
      dockerfile: Dockerfile
    volumes:
      - ./src:/src
    environment:
      - TZ=Asia/Tokyo
    restart: always

フォルダ構成は以下の通りとします。

Docker起動後、srcフォルダ内で作業することを想定しています。

.
├── docker
│   └── Dockerfile
├── docker-compose.yml
└── src

あとは、以下のようにして起動し、作業を開始します。

docker compose up -d --build
docker compose exec now_cast_amedas bash # コンテナに入る

cd /src # 作業ディレクトリに移動

サンプルデータ

気象庁が提供しているサンプルデータを利用します。

www.data.jma.go.jp
www.data.jma.go.jp

この中の 高解像度降水ナウキャスト のサンプルをダウンロードします。

ダウンロード後、srcフォルダに配置します。

これでサンプルデータの準備は完了です。

中身の確認

zipファイルの解凍

まず、zipファイルを解凍します。

unzip highnowc.zip 

Archive:  highnowc.zip
  inflating: Z__C_RJTD_20210817090000_NOWC_GPV_Ggis0p25km_Pri60lv_Aper5min_FH0000-0030_grib2.bin.gz  
  inflating: Z__C_RJTD_20210817090000_NOWC_GPV_Ggis0p25km_Prr05lv_Aper5min_FH0000-0030_grib2.bin.gz 

2つのファイルが生成されます。

それぞれの違いは以下に記載されています。

www.data.jma.go.jp
www.data.jma.go.jp
ファイル名説明
Z__C_RJTD_xx_NOWC_GPV_Ggis0p25km_Pri60lv_Aper5min_FH0000-0030_grib2.bin.gz瞬間的な降水強度である「高解像度降水ナウキャスト」
Z__C_RJTD_xx_NOWC_GPV_Ggis0p25km_Prr05lv_Aper5min_FH0000-0030_grib2.bin.gz5分間降水量である「高解像度降水ナウキャスト(5分間)」

気象庁が公開している雨雲サイト と同じ雨雲データにしたい」場合は、Z__C_RJTD_xx_NOWC_GPV_Ggis0p25km_Pri60lv_Aper5min_FH0000-0030_grib2.bin.gz になると思われます。

降水強度を出しているとのことなので、おそらくこれで合っているでしょう。

気象庁の雨雲サイト

今回は、降水強度のデータを利用します。

bin.gzファイルの解凍

以下のコマンドでファイルを解凍します。

gunzip Z__C_RJTD_20210817090000_NOWC_GPV_Ggis0p25km_Pri60lv_Aper5min_FH0000-0030_grib2.bin.gz 

解凍すると、Z__C_RJTD_20210817090000_NOWC_GPV_Ggis0p25km_Pri60lv_Aper5min_FH0000-0030_grib2.binというファイルが生成されます。

中身を見てみる

wgrib2コマンドを使って中身を見てみます。

多くの文字が出力されます。

この部分はフォーマットが決まっているようですが、何がどの値に対応しているのかが分かりにくいです。
参考:フォーマットは 仕様書のp4の別紙1 に記載されています。

wgrib2 Z__C_RJTD_20210817090000_NOWC_GPV_Ggis0p25km_Pri60lv_Aper5min_FH0000-0030_grib2.bin

1.1:0:d=2021081709:var discipline=0 center=34 local_table=1 parmcat=1 parm=203:surface:-5-0 min ??? fcst:
1.2:0:d=2021081709:var discipline=0 center=34 local_table=1 parmcat=1 parm=203:surface:-5-0 min ??? fcst:
1.3:0:d=2021081709:var discipline=0 center=34 local_table=1 parmcat=1 parm=203:surface:-5-0 min ??? fcst:

~~ 略

1.24819:0:d=2021081709:var discipline=0 center=34 local_table=1 parmcat=1 parm=203:surface:55-60 min ??? fcst:
1.24820:0:d=2021081709:var discipline=0 center=34 local_table=1 parmcat=1 parm=203:surface:55-60 min ??? fcst:
1.24821:0:d=2021081709:var discipline=0 center=34 local_table=1 parmcat=1 parm=203:surface:55-60 min ??? fcst:
1.24822:0:d=2021081709:var discipline=0 center=34 local_table=1 parmcat=1 parm=203:surface:55-60 min ??? fcst:

ChatGPTによると、以下の通りです。

確かに仕様書を見ると、このような記述があります。

surfaceの値を見る限り、5分刻みで60分までのデータが含まれているようです。

各フィールドの意味はこうなります:

1.1
→ これは wgrib2 が付けたメッセージ番号(インデックス)。-inv でインデックスファイルを作ると、この番号で指定して取り出せます。

0
→ GRIB2 メッセージ内でのバイトオフセット。

d=2021081709
→ データ時刻(基準時刻)= 2021年8月17日 09:00 UTC。

var discipline=0
→ discipline=0 は「気象」分野(GRIB2 discipline 0)。

center=34
→ 気象庁 (Tokyo, Japan Meteorological Agency) のセンターID。

local_table=1 parmcat=1 parm=203
→ ここがポイント。

parmcat=1 → 降水関連(precipitation)。

parm=203 → JMA 独自のパラメータ番号。
→ つまり「JMA独自の降水量関連データ(ナウキャスト用変数)」です。

surface
→ 高度座標(ここでは地表)。

-5-0 min ??? fcst
→ 予報時間範囲。

-5-0 min は「直近5分間の積算値」などを意味している可能性が高いです。

??? fcst は wgrib2 が変数定義を知らないためにラベルを付けられず「???」と出している状態。

ちなみに、ここで出力された内容には、降水強度の具体的な数値は含まれていません。あくまで概要のみです。

仕様書のp4の別紙1 に従って変換することで、具体的な数値が分かります。

数値を変換する方法としては、バイナリの中身を直接解析するか、netCDF形式に変換するかのどちらかになるでしょう。

前者の方法は以下のサイトを参考にしてください。

GRIB2からマップタイルを出力したい - Qiita
qiita.com
image

今回は後者のnetCDFに変換する方法を採用します。

netCDFへの変換

wgrib2コマンドでnetCDF形式に変換します。

ちなみに、netCDFとは、気象、海洋、地球科学分野でよく用いられる科学データファイル形式の一つです。
参考: https://ja.wikipedia.org/wiki/NetCDF

まず、1つのIDずつ(専門用語では1メッセージと呼ぶようです)指定して変換を行います。
注意:まとめて一気に変換することはできません。

補足: まとめて一気に変換できない理由

実際に試してみると、以下のように、1つのグリッドタイプしか変換できないというエラーが発生します。

wgrib2 nowc_grib2_-5-0.bin -netcdf out.nc # あらかじめ -5-0 min のデータを抽出したファイルを指定

1:0:d=2021081709:var discipline=0 center=34 local_table=1 parmcat=1 parm=203:surface:-5-0 min ??? fcst:

*** FATAL ERROR: netcdf: only one grid type at a time ***

*** arg list to wgrib2(..): wgrib2 nowc_grib2_-5-0.bin -netcdf out.nc

とりあえず、先頭のデータをnetCDFに変換してみます。

wgrib2 Z__C_RJTD_20210817090000_NOWC_GPV_Ggis0p25km_Pri60lv_Aper5min_FH0000-0030_grib2.bin -d 1.1 -netcdf cdf_1.nc

1.1:0:d=2021081709:var discipline=0 center=34 local_table=1 parmcat=1 parm=203:surface:-5-0 min ??? fcst:

このファイルをncdumpコマンドで中身を見てみます。

ncdump -h cdf_1.nc 

netcdf cdf_1 {
dimensions:
 latitude = 40 ;
 longitude = 40 ;
 time = UNLIMITED ; // (1 currently)
variables:
 double latitude(latitude) ;
  latitude:units = "degrees_north" ;
  latitude:long_name = "latitude" ;
 double longitude(longitude) ;
  longitude:units = "degrees_east" ;
  longitude:long_name = "longitude" ;
 double time(time) ;
  time:units = "seconds since 1970-01-01 00:00:00.0 0:00" ;
  time:long_name = "verification time generated by wgrib2 function verftime()" ;
  time:reference_time = 1629190800. ;
  time:reference_time_type = 1 ;
  time:reference_date = "2021.08.17 09:00:00 UTC" ;
  time:reference_time_description = "analyses, reference date is fixed" ;
  time:time_step_setting = "auto" ;
  time:time_step = 0. ;
 float var0_1_203_surface(time, latitude, longitude) ;
  var0_1_203_surface:_FillValue = 9.999e+20f ;
  var0_1_203_surface:short_name = "var0_1_203_surface" ;
  var0_1_203_surface:long_name = "desc" ;
  var0_1_203_surface:level = "surface" ;
  var0_1_203_surface:units = "unit" ;

// global attributes:
  :Conventions = "COARDS" ;
  :History = "created by wgrib2" ;
  :GRIB2_grid_template = 0 ;
}

正直、まだよく分かりません。

とりあえず、画像化してみます。

tiff形式に変換してみる

netCDFファイルからGeoTIFF(地理情報を含む画像ファイル)に変換してみます。

gdal_translate -of GTiff cdf_1.nc out_1.tif

Input file size is 40, 40
0...10...20...30...40...50...60...70...80...90...100 - done.

Macのプレビューで開いてみます。

GeoTiff変換時にInput file size is 40, 40と表示された通り、40x40ピクセルの画像が生成されています。

画像画像情報

正直なところ、黒い画像が生成されるだけで、何を示しているのか判別できません。

ちなみに、他のIDについても同様に変換すると、何らかの模様が現れます。

色付けしてみる

白黒のtiff画像を色付けしてみます。

あらかじめ色付けのルールを定義しておきます。

テキストファイル(color.txt)に以下のように記述します。

データの欠損値(nv9.999e+20)は透過とし、それ以外は降水強度に応じて色を変更します。

nv      0 0 0 0
9.9990003e+20   0 0 0 0
9.999e+20       0 0 0 0
-9999   0 0 0 0
0       255 255 255 0
1       0   128 255 255
5       0   255 255 255
10      0   255 0   255
20      255 255 0   255
30      255 128 0   255
50      255 0   0   255
80      128 0   0   255

それでは、色付けを実行します。

gdaldem color-relief out_1.tif color.txt out_color_1.tiff -alpha

0...10...20...30...40...50...60...70...80...90...100 - done.

色付けしたtiffファイルを開いてみます。

白い四角形が生成されます。

色付け結果

gdalinfoコマンドで情報を確認してみます。

座標情報が含まれており、これによって画像を地図上に配置できるようです。

gdalinfo out_color_1.tiff 

Driver: GTiff/GeoTIFF
Files: out_color_1.tiff
Size is 40, 40
Origin = (140.000000000000000,46.666666666666671)
Pixel Size = (0.012499999999989,-0.008333333333333)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 140.0000000,  46.6666667) 
Lower Left  ( 140.0000000,  46.3333333) 
Upper Right ( 140.5000000,  46.6666667) 
Lower Right ( 140.5000000,  46.3333333) 
Center      ( 140.2500000,  46.5000000) 
Band 1 Block=40x40 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=40x40 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=40x40 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=40x40 Type=Byte, ColorInterp=Alpha

他のIDについても同様に色付けすると、適切に色がつけられていることが分かります。

色付け結果

pmtiles形式に変換してみる

地図上に配置するためにpmtiles形式に変換してみます。

pmtilesに変換すると、タイル画像を1つのファイルにまとめることができます。

実務上では、すべてのタイルを結合して1つのpmtilesを作成することになると思いますが、まずは1つだけで試してみます。

まず、gdal_edit.pyコマンドでtiffファイルに座標参照系を設定します。

EPSG:4326を設定します。
※これを設定しないと、地図上に画像を配置できません。

 gdal_edit.py -a_srs EPSG:4326 out_color_1.tiff

次に、rio pmtilesコマンドでpmtilesに変換します。

webp形式で圧縮します。

rio pmtiles out_color_1.tiff out_color_1.pmtiles --format WEBP --resampling bilinear

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 7/7 [00:00<00:00, 166.54it/s]

無事に生成されたので、pmtiles viewerで開いてみます。

何も見えません…。下記画像の赤枠にあるとのことです(透過画像なので見えないのだと思われます)。

pmtiles viewer の結果1

別のGeoTIFF(メッセージID: 1354)で試してみます。

適切に表示されました。

pmtiles viewer の結果2

ここまでできれば、あとは複数のGeoTIFFを結合してpmtilesにすれば、地図上に配置できます。

ある時間帯のGeoTIFFを結合してpmtilesにする

ここまで長々と説明してきましたが、本題である「ある時間帯のGeoTIFFを結合してpmtilesにする」を試してみます。

grib2の実行結果から分かるように、データは5分刻みで60分先まで含まれています。

  • -5-0 min
  • 0-5 min
  • 5-10 min
  • 10-15 min
  • 15-20 min
  • 20-25 min
  • 25-30 min
  • 30-35 min
  • 35-40 min
  • 40-45 min
  • 45-50 min
  • 50-55 min
  • 55-60 min

今回は-5-0 minのデータを結合してpmtilesにしてみます。

ソースコードやフォルダ構成は折りたたんでおきます。

フォルダ構成/ソースコード

フォルダ構成

ナウキャストのサンプルデータをsrc/nowcastフォルダに配置します。

そして、pmtilesを生成するスクリプトをsrc/main.pyに置きます。

.
├── docker
│   └── Dockerfile
├── docker-compose.yml
└── src
    ├── color.txt
    ├── main.py
    └── nowcast
        └── highnowc.zip

ソースコード

行っている処理は以下の通りです。

  1. -5-0 minの範囲にあるメッセージを1件ずつnetCDF形式に変換
  2. netCDFファイルをGeoTIFF形式に変換
  3. すべてのGeoTIFFファイルを結合し、VRT(Virtual Raster)ファイルを作成
  4. VRTファイルをGeoTIFF形式に変換
  5. GeoTIFFファイルに色付けを行う
  6. 色付けしたGeoTIFFファイルをpmtiles形式に変換

手順1は、3,546件あるため、並列処理で高速化しています。

import os
import zipfile
import glob
import subprocess
import gzip
import shutil
from netCDF4 import Dataset
import numpy as np
from concurrent.futures import ProcessPoolExecutor, as_completed

ZIP_PATH = "./src/nowcast/highnowc.zip"

GRIB2_FILE_PATTERN = "Z__C_RJTD_*_NOWC_GPV_Ggis0p25km_Pri60lv_Aper5min_FH0000-0030_grib2"
GRIB2_GZIP_FILE = f"{GRIB2_FILE_PATTERN}.bin.gz"
GRIB2_FILE = "nowc_grib2.bin"

TMP_DIR = f"./src/tmp/{os.urandom(8).hex()}"
GRIB2_DIR = os.path.join(TMP_DIR, "grib2")
TILE_DIR = os.path.join(TMP_DIR, "tiles")

OUTPUT_DIR = "./src/output"
COLOR_TXT = "./src/color.txt"

# CPU コア数を取得
NUM_WORKERS = os.cpu_count()

def run_command(cmd, **kwargs):
    subprocess.run(cmd, check=True, **kwargs)


def extract_zip(zip_path, extract_to):
    """zip の解凍をする関数"""
    with zipfile.ZipFile(zip_path, "r") as zip_ref:
        zip_ref.extractall(extract_to)
    print(f"Extracted {zip_path} to {extract_to}")
    return extract_to


def extract_gzip(gzip_path, extract_to):
    """gzip の解凍をする関数"""
    with gzip.open(gzip_path, "rb") as f_in:
        with open(extract_to, "wb") as f_out:
            shutil.copyfileobj(f_in, f_out)
    print(f"Extracted {gzip_path} to {extract_to}")
    return extract_to


def get_grib2_messages(grib_file, target_period):
    """grib2 から対象時刻を取得する関数"""
    output_file = f"{grib_file[:-4]}_{target_period}.bin"
    run_command(
        [
            "wgrib2",
            grib_file,
            "-match",
            f":surface:{target_period} min",
            "-grib_out",
            output_file,
        ]
    )
    return output_file


def process_message(msg_info):
    """メッセージを処理する関数"""
    i, msg_num, grib_file, tile_dir, nodata_value = msg_info
    nc_tmp = f"{tile_dir}/rain_{i+1}.nc"
    tif_tmp = f"{tile_dir}/rain_{i+1}.tif"
    try:
        run_command(
            ["wgrib2", grib_file, "-d", msg_num.split(":")[0], "-netcdf", nc_tmp]
        )
        with Dataset(nc_tmp) as ds:
            if "var0_1_214_surface" in ds.variables:
                rain_data = ds.variables["var0_1_214_surface"].data
                rain_data = np.nan_to_num(rain_data, nan=nodata_value)
                ds.variables["var0_1_214_surface"].data = rain_data

            run_command(
                [
                    "gdal_translate",
                    "-of",
                    "GTiff",
                    "-a_nodata",
                    str(nodata_value),
                    nc_tmp,
                    tif_tmp,
                ]
            )
        os.remove(nc_tmp)
        return (i, None)
    except Exception as e:
        print(f"Error processing {nc_tmp}: {e}")
        try:
            os.remove(nc_tmp)
        except Exception:
            pass
        return (i, e)


def create_grib2_to_tiffs(grib_file):
    """grib2 から NetCDF ファイルを作成する関数"""
    msg_list = (
        subprocess.check_output(["wgrib2", grib_file, "-s"]).decode().splitlines()
    )
    nodata_value = 9.999e20
    tile_dir = TILE_DIR

    msg_infos = [
        (i, msg_num, grib_file, tile_dir, nodata_value)
        for i, msg_num in enumerate(msg_list)
    ]

    with ProcessPoolExecutor(max_workers=NUM_WORKERS) as executor:
        futures = [executor.submit(process_message, msg_info) for msg_info in msg_infos]
        for future in as_completed(futures):
            i, err = future.result()
            if err:
                print(f"Failed to process message #{i}: {err}")


def create_pmtiles_from_tiffs(tif_dir, output_dir):
    """tiff ファイルを結合して pmtiles に変換する関数"""
    vrt_file = os.path.join(output_dir, "tiles.vrt")
    tif_list = glob.glob(os.path.join(tif_dir, "*.tif"))
    run_command(["gdalbuildvrt", vrt_file] + tif_list)

    tif_file = os.path.join(output_dir, "tiles.tif")
    run_command(["gdal_translate", "-of", "GTiff", vrt_file, tif_file])

    tif_colored_file = os.path.join(output_dir, "tiles_colored.tif")
    run_command(
        [
            "gdaldem",
            "color-relief",
            tif_file,
            COLOR_TXT,
            tif_colored_file,
            "-alpha",
        ]
    )

    run_command(["gdal_edit.py", "-a_srs", "EPSG:4326", tif_colored_file])

    pmtiles_file = os.path.join(output_dir, "tiles.pmtiles")
    run_command(
        [
            "rio",
            "pmtiles",
            tif_colored_file,
            pmtiles_file,
            "--format",
            "WEBP",
            "--resampling",
            "bilinear",
        ]
    )


def main():
    period = "-5-0"

    # tmp フォルダの作成
    os.makedirs(TMP_DIR, exist_ok=True)
    os.makedirs(GRIB2_DIR, exist_ok=True)
    os.makedirs(TILE_DIR, exist_ok=True)
    os.makedirs(os.path.join(OUTPUT_DIR, period), exist_ok=True)

    # zip 解凍
    extract_zip(ZIP_PATH, GRIB2_DIR)

    # gzip 解凍
    target_gz_files = glob.glob(os.path.join(GRIB2_DIR, GRIB2_GZIP_FILE))
    if not target_gz_files:
        print(f"No gzip file matching {GRIB2_GZIP_FILE} found in {GRIB2_DIR}")
        return
    target_gz_file = target_gz_files[0]

    grib_file = extract_gzip(target_gz_file, os.path.join(GRIB2_DIR, GRIB2_FILE))

    # 対象時刻の GRIB2 ファイルを取得
    target_period_grib2_file = get_grib2_messages(grib_file, period)

    # grib2 から tiff ファイルを作成
    create_grib2_to_tiffs(target_period_grib2_file)

    # 複数 tiff を 結合して pmtiles に変換
    tif_files = glob.glob(os.path.join(TILE_DIR, "*.tif"))
    if not tif_files:
        print("No TIFF files found to convert to PMTiles.")
        return
    create_pmtiles_from_tiffs(TILE_DIR, os.path.join(OUTPUT_DIR, period))

if __name__ == "__main__":
    main()

起動方法は以下の通りです。

cd /
python src/main.py

Extracted ./src/nowcast/highnowc.zip to ./src/tmp/ce5f9fefbcc3d188/grib2
Extracted ./src/tmp/ce5f9fefbcc3d188/grib2/Z__C_RJTD_20210817090000_NOWC_GPV_Ggis0p25km_Pri60lv_Aper5min_FH0000-0030_grib2.bin.gz to ./src/tmp/ce5f9fefbcc3d188/grib2/nowc_grib2.bin
~~~ 略 ~~~

ソースコードを実行すると、ログに以下のようなエラーが出てきています。

このエラーはデータの欠損が原因で発生しています。無視しても問題ないです。

Failed to process message #1814: NetCDF: Attribute not found

Error processing ./src/tmp/ce5f9fefbcc3d188/tiles/rain_1817.nc: NetCDF: Attribute not found

コードを実行すると、 ./src/output/-5-0/ に生成物ができています。

  • tiles_colored.tif ・・・ 色付けした GeoTiff
  • tiles.pmtiles ・・・ pmtiles ファイル
  • tiles.tif ・・・ 色付けしていない 結合した GeoTiff。
  • tiles.vrt ・・・ 結合した GeoTiff の情報
色付けした GeoTiff (tiles_colored.tif)
pmtiles viewer で tiles.pmtiles を開いた結果
いい感じにできているっぽい

まとめ

高解像度降水ナウキャストのサンプルデータを使って、ある時間帯の GeoTiff を結合して pmtiles にする方法を紹介しました。

どなたかの参考になれば幸いです。

おかしなところあれば、X(旧 Twitter)で直接教えてください。

x.com
x.com

github

この記事では細かく書きましたが、とりあえず動くコードがほしいという人は github に置いておきました。

GitHub - Momijinn/SampleNowCastHighnowc: 気象庁が提供している高解像度降水ナウキャストを画像化するサンプル
github.com
image