GCOM-Cのデータで全球画像を作ってみる (その2)

その1から2年近く立ってしまいましたが、GCOM-Cの陸域反射率のファイルを使用して全球タイルマップを作成しました。 COGファイルをソースとしたタイルサーバを作成しました。

全球GeoTIFFファイルをつくる
今回作成したアプリのフロー
x,y,zからタイル画像の緯度経度の範囲を導出
STAC catalogから緯度経度の範囲に該当するCOG検索し、そのリストを取得する
アプリ実行
参考文献

全球GeoTIFFファイルをつくる

以前紹介した方法でGCOM-C L2 統計陸域反射率のファイルをG-Portalからダウンロードし、COGファイルを作成します。そのCOGファイルをS3に保存します。

GCOM-C L2のHDF5 データを地図投影までしてみる。(その2) - たぶん動く...

このままCOGのリストを以前紹介した方法にリストとして渡してもタイルサーバは動きますが、rio-tilerが画像出力する場所に相当するファイルを検索するのに時間がかかり、レスポンスが非常に長くなり使い物になりません。
そこでSTACのリスト (STAC Catalog)を作成し、そのリストからファイルを検索するようにします。

今回作成したアプリのフロー

アクティビティ図を作成してみました。このうち、①、②の部分について解説します。他は以前紹介した方法とほとんど変わりません。

x,y,zからタイル画像の緯度経度の範囲を導出

ブラウザやQGISからアクセスすると、Tile Map Serviceで使われる地図タイルのリクエストがサーバに届きます。リクエストされるタイルからその画像の範囲の緯度経度を計算します。

計算方法のいい文献がなかったので、ChatGPTに聞きました。計算式は以下です。タイルマップでxyzが与えられたとき、左上の緯度経度(lat,lon)は以下で求まります。

$n = {2}^{z}$

$lon = ( \dfrac{x}{n} \times 360.0) - 180.0$

$lat = \arctan(\sinh( \pi \times (1 - 2 \times \dfrac{y}{n})))$

これを応用してタイルの右下の緯度経度 (lat, lon)は以下で求まります。

$n = {2}^{z}$

$lon = ( \dfrac{x+1}{n} \times 360.0) - 180.0$

$lat = \arctan(\sinh( \pi \times (1 - 2 \times \dfrac{y+1}{n})))$

プログラムで計算する際は、degreeとradianの扱いに注意してください。

STAC catalogから緯度経度の範囲に該当するCOG検索し、そのリストを取得する

今回制作したアプリのミソになります。 COGからSTACとそのSTAC catalogを作成しました。STAC catalogからタイル画像を作成するために必要なCOGを取得します。 COGを複数枚ホストする方法と同じように、全球分のCOGファイルをassetsに使用してもいいのですが、タイル画像作成に数分かかります。今回の方法はより高速にタイル画像を作成できます。

pystacのインストール

pystacをインストールすると、stac関連をコマンドラインで処理することができます。

pip install pystac

COGからstacのメタデータを作成する

stac create-itemでs3に配置したCOGへのパスを指定すると、stac形式でメタデータが出力されます。

stac create-item s3://COG-PATH/COG-FILE

{"type": "Feature", "stac_version": "1.0.0", "id": "GC1SG1_20230407D08D_T0112_L2SG_RGB358", "properties": {"proj:wkt2": "PROJCS[\"Sphere_Sinusoidal\",GEOGCS[\"Unknown datum based upon the Authalic Sphere\",DATUM[\"Not_specified_based_on_Authalic_Sphere\",SPHEROID[\"Sphere\",6371000,0,AUTHORITY[\"EPSG\",\"7035\"]],AUTHORITY[\"EPSG\",\"6035\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Sinusoidal\"],PARAMETER[\"longitude_of_center\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]", "proj:transform": [231.65635828469235, 0.0, -6671703.118599138, 0.0, -231.65635828469235, 8895604.158132184], "proj:shape": [4800, 4800], "datetime": "2023-05-22T13:19:39.089577Z"}, "geometry": {"type": "Polygon", "coordinates": [[[-146.190938, 70.000079], [72.058581, 80.00009], [14.470298, 80.00009], [-175.429125, 70.000079], [-146.190938, 70.000079]]]}, "links": [], "assets": {"data": {"href": "s3://COG-PATH/COG-FILE", "roles": ["data"]}}, "bbox": [-175.429125, 70.000079, 72.058581, 80.00009], "stac_extensions": ["https://stac-extensions.github.io/projection/v1.1.0/schema.json"]

これをCOGファイル1つずつに対し作成し、S3へ保存します。作成した全部のCOGのリストを作成します。これがSTACのカタログに相当します。catalog.jsonというファイル名で作成しました。これもs3に保存します。 STACコマンドで作れるはずですが、よくわからなかったので手打ちで作成しました。rootにはcatalog.jsonを指定します。

{
    "id": "sample-catalog",
    "stac_version": "1.0.0",
    "description": "Simple STAC catalog.",
    "links": [
        {
            "rel": "root",
            "href": "s3://path-to-s3-bucket/SGLI/L2/catalog.json",
            "type": "application/json"
        },
        {
            "rel": "item",
            "href": "s3://path-to-s3-bucket/SGLI/L2/SGLI/GC1SG1_20230322D08D_T0320_L2SG_RGB468.json",
            "type": "application/json"
...
        {
            "rel": "item",
            "href": "s3://path-to-s3-bucket/SGLI/L2/SGLI/GC1SG1_20230322D08D_T0829_L2SG_RGB468.json",
            "type": "application/json"
        }
    ]
}

STAC catalogから緯度経度の範囲に該当するCOG検索し、リストを取得する

pystac-clientモジュールを使用します。

from pystac_client import Client

catalog= Client.open('https://path-to-s3-bucket/catalog.json')  <- STAC catalogを呼び出す
child_link=catalog.get_links('item') <- COG個々のSTACのリンクをリストで持つ

 aoi=bbox_from_xyz(x,y,z)    <- 先述の方法でxyzからaoiを導出
    item_links=[]
    for item in range(len(child_link)):
        item_bbox=child_link[item].resolve_stac_object(catalog).target.bbox
        cog_bbox=pyproj.aoi.BBox(item_bbox[0],item_bbox[1],item_bbox[2],item_bbox[3])   <- aoiがstacの持っているCOGファイルのaoiと重なるか検索する
        if aoi.intersects(cog_bbox) :  <- aoiが重なったらリストに追加する
            stac_href=child_link[item].resolve_stac_object(catalog).href
            item_links.append(child_link[item].resolve_stac_object(catalog).target.assets["data"].href)

アプリ実行

画像を扱うので、ある程度のコンピューティングスペックが必要です。AWSのEC2で実行する場合、t3a.microでは途中でハングしてしまいます。 t3a.medium以上のサイズをお勧めします。

uvicorn app:app --host 0.0.0.0 --timeout-keep-alive 100 --workers 8 --limit-max-requests 2000

実行ログを見てみます。 REQUESTがサーバがリクエストを受け取った時刻、send pngがタイル画像を送信した時刻になります。 zoomレベル9では、おおよそ5秒くらいで返せています。
いろいろいじってみたのですが、表示範囲が大きいと画像生成に時間がかかり、レスポンスが悪くなるようです。
次作成するときは、STACサーバをちゃんと立ててアプリを作成してみたいです。

REQUEST     9 447 202 2023-05-22 14:02:52.978526
REQUEST     9 446 202 2023-05-22 14:02:52.989101
REQUEST     9 447 203 2023-05-22 14:02:53.000071
REQUEST     9 447 204 2023-05-22 14:02:53.007933
REQUEST     9 446 203 2023-05-22 14:02:53.038547
REQUEST     9 447 201 2023-05-22 14:02:53.042972
send png    9 447 203 2023-05-22 14:02:57.003382
INFO:     127.0.0.1:47776 - "GET /9/447/203.png HTTP/1.1" 200 OK
REQUEST     9 446 201 2023-05-22 14:02:57.050804
send png    9 447 202 2023-05-22 14:02:57.603700
INFO:     127.0.0.1:47758 - "GET /9/447/202.png HTTP/1.1" 200 OK
REQUEST     9 446 204 2023-05-22 14:02:57.661303
send png    9 446 202 2023-05-22 14:02:58.882789
INFO:     127.0.0.1:47768 - "GET /9/446/202.png HTTP/1.1" 200 OK
REQUEST     9 445 202 2023-05-22 14:02:58.929019

参考文献

rio-tiler
PySTAC Documentation — pystac 1.7.3 documentation