雲だけのラスタータイルが欲しいなと思いました。なるべく地球全体のものを。
以下のデータがありました。
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'
次に、雲の部分だけ欲しいのでこのデータの黒い部分を透過させます。これは以下のコマンドで対応しました。
# 各ピクセルの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'
このような見た目になります。
これを、雲のない衛星写真のレイヤーと組み合わせることで、地表レイヤーと雲レイヤーに分けて表示できます
衛星写真のレイヤーはUSGSImageryOnlyを使用してます。
これを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,
});
こんな感じになります。
以下のようにズーム式を使ってraster-opacityを可変するようにします。
{
id: 'cloud_layer',
source: 'cloud',
type: 'raster',
paint: {
'raster-opacity': [
'interpolate',
['linear'],
['zoom'],
4,
1.0,
5,
0,
],
},
},
すると、ズームレベルによって雲の表示が切り替わる表現ができます。これはGoogleMapsの地球表示のときにやっている表現を参考にしました。
DEMOはこちらになります。