環境:Houdini 20.5.487

Mapboxから生成した立体地図を球面に投影する。
処理の流れ

Mapboxから生成した衛星画像から地形を用意する。
Mapboxのハイトマップから地形を生成する


地球は半径約6,371kmの球とする。緯度と経度からポイントを球の座標へ変換する。中心座標のフォワードベクトルとアップベクトルを計算し行列をつくり、逆行列でポイントを原点へ戻す。

回転による座標変換の図と計算式。

球状に曲がった状態で原点に戻る。
コード
//
// 平面に配置されている座標を地球の曲面に合わせて変換させる
// RunOver: Points
//
// ハイトマップ四方の長さ
float width = `chs("../Coordinates/width")`;
float height = `chs("../Coordinates/height")`;
// 四方の緯度と経度
float west = `chs("../Coordinates/west")`;
float east = `chs("../Coordinates/east")`;
float north = `chs("../Coordinates/north")`;
float south = `chs("../Coordinates/south")`;
// 地球の半径
float radiusEarth = 6371000;
// 緯度・経度・高度を調べる
// 標高
float h = @P.y;
// 緯度
float center_lat = (north - south) / 2 + south;
float lat = @P.z * -1 / (height/2) * (north - south) / 2 + center_lat;
// 経度
float center_lon = (east - west) / 2 + west;
float lon = @P.x / (width/2) * (east - west) / 2 + center_lon;
// 球面上に座標変換する
float y = sin(radians(lat));
float x = cos(radians(lat)) * cos(radians(lon)*-1);
float z = cos(radians(lat)) * sin(radians(lon)*-1);
vector pos = set(x, y, z) * (radiusEarth + h);
// 中心点の座標
y = sin(radians(center_lat));
x = cos(radians(center_lat)) * cos(radians(center_lon)*-1);
z = cos(radians(center_lat)) * sin(radians(center_lon)*-1);
// 中心点の行列をつくる
vector center = set(x, y, z) * radiusEarth;
vector N = normalize(center);
vector up = set(0, 1, 0);
vector left = cross(up, N);
up = normalize(cross(N, left));
matrix world = maketransform(N, up, center);
matrix inverse = invert(world);
// 行列を90°回転させる
float angle = PI/2*-1;
vector axis = set(1, 0, 0);
rotate(inverse, angle, axis);
pos *= inverse;
@P = pos;