土曜日, 7月 5, 2025
土曜日, 7月 5, 2025
- Advertisment -
ホームニューステックニュース雲だけのラスタータイルを作成する #GIS - Qiita

雲だけのラスタータイルを作成する #GIS – Qiita



雲だけのラスタータイルを作成する #GIS - Qiita

雲だけのラスタータイルが欲しいなと思いました。なるべく地球全体のものを。
以下のデータがありました。

image.png

NASA/Goddard Space Flight Center, Image by Reto Stöckli

このデータセットは、NASAの「Blue Marble」データセットの雲画像部分です。地球全体の雲システムを継ぎ目なく表現したデータで、白い部分が雲の密度を表しています。ただし、提供されているTIFF画像にはGIS座標情報が含まれておらず、通常のRGB写真画像として配布されています。

これをgdalを使って位置情報データに変換して可視化します。

以下のコマンドでEPSG:4326の座標を付与します。これでQGISにデータが乗るようになります。

gdal_translate \
    -a_srs EPSG:4326 \
    -a_ullr -180 90 180 -90 \
    'cloud_combined_8192.tif' \
    'cloud_combined_8192_georef.tif'

image.png

次に、雲の部分だけ欲しいのでこのデータの黒い部分を透過させます。これは以下のコマンドで対応しました。

# 各ピクセルのR・G・B値の最大値を計算し、透明度マスクを作成
# 白い雲部分は不透明、黒い背景部分は透明になる
gdal_calc.py -A '/Users/satoshi/Downloads/cloud_combined_8192_georef.tif' --A_band=1 \
  -B 'cloud_combined_8192_georef.tif' --B_band=2 \
  -C 'cloud_combined_8192_georef.tif' --C_band=3 \
  --outfile='clouds_transparent.tif' \
  --calc="maximum(maximum(A, B), C)" \
  --creation-option="COMPRESS=LZW" \
  --overwrite

# 元のRGB画像とアルファチャンネルを結合してRGBA画像を作成
gdal_merge.py -separate -o 'clouds_rgba.tif' \
  'cloud_combined_8192_georef.tif' \
  'clouds_transparent.tif'

このような見た目になります。

image.png

これを、雲のない衛星写真のレイヤーと組み合わせることで、地表レイヤーと雲レイヤーに分けて表示できます
衛星写真のレイヤーはUSGSImageryOnlyを使用してます。

Kapture 2025-07-04 at 14.30.18.gif

これをWebGISで表示するためにラスタータイル(pmtiles)に変換します。こちらの記事を参考にしました。

pmtilesに変換(go-pmtilesの使用)

#EPSG:3857に変換
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 clouds_rgba.tif clouds_rgba_mercator.tif

#mbtilesに変換
gdal_translate -of MBTILES clouds_rgba_mercator.tif cloud.mbtiles -r bilinear -co BLOCKSIZE=512 -co TILE_FORMAT=WEBP

#低いズームレベルのタイルを生成する
gdaladdo -r bilinear cloud.mbtiles

#go-pmtilesのバイナリファイルを使用してmbtilesをpmtileに変換
./pmtiles convert cloud.mbtiles cloud.pmtiles

MapLIbre GL JSを使って表示してみます。

import maplibregl from 'maplibre-gl';
import 'maplibre-gl/dist/maplibre-gl.css';

import { Protocol } from 'pmtiles';

// PMTilesのプロトコルを設定
const pmtilesProtocol = new Protocol();
maplibregl.addProtocol('pmtiles', pmtilesProtocol.tile);

// 地図の表示
const map = new maplibregl.Map({
    container: 'map',
    style: {
        projection: {
            type: 'globe',
        },
        version: 8,
        glyphs: 'https://demotiles.maplibre.org/font/{fontstack}/{range}.pbf',
        sources: {
            // USGSの衛星画像タイル
            USGSImageryOnly: {
                type: 'raster',
                tiles: ['https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}'], // タイルのURL
                tileSize: 256,
                maxzoom: 10,
                attribution: 'USGS',
            },
            // 雲レイヤーのPMTiles
            cloud: {
                type: 'raster',
                url: `pmtiles://./cloud.pmtiles`,
                attribution: 'NASA/Goddard Space Flight Center, Image by Reto Stöckli',
            },
        },
        layers: [
            {
                id: 'USGSImageryOnly_layer',
                source: 'USGSImageryOnly',
                type: 'raster',
            },
            {
                id: 'cloud_layer',
                source: 'cloud',
                type: 'raster',
                paint: {
                    'raster-opacity': 1.0,
                },
            },
        ],
        sky: {
            'sky-color': '#199EF3',
            'sky-horizon-blend': 0.5,
            'horizon-color': '#ffffff',
            'horizon-fog-blend': 0.5,
            'fog-color': '#0000ff',
            'fog-ground-blend': 0.5,
            'atmosphere-blend': ['interpolate', ['linear'], ['zoom'], 0, 1, 10, 1, 12, 0],
        },
    },

    center: [139.477, 35.681],
    zoom: 2,
});

こんな感じになります。

image.png

以下のようにズーム式を使ってraster-opacityを可変するようにします。

 {
    id: 'cloud_layer',
    source: 'cloud',
    type: 'raster',
    paint: {
        'raster-opacity': [
            'interpolate',
            ['linear'],
            ['zoom'],
            4,
            1.0,
            5,
            0,
        ],
    },
},

すると、ズームレベルによって雲の表示が切り替わる表現ができます。これはGoogleMapsの地球表示のときにやっている表現を参考にしました。

Kapture 2025-07-04 at 14.42.06.gif

DEMOはこちらになります。





Source link

Views: 0

RELATED ARTICLES

返事を書く

あなたのコメントを入力してください。
ここにあなたの名前を入力してください

- Advertisment -