地図の球面投影

環境: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;
タイトルとURLをコピーしました