日々適当

hibitekitou

正距方位図法に変換

xsi |2015-07-07
地図の中心からの各地点への距離と方向が正しくなるという図法で、日本地図[国土地理院]だと北緯35度39分29秒、東経139度44分29秒という位置が中心(主点)となるようです。

CG上で地球のテクスチャとしてよく使われる正距円筒図法は、経度、緯度をそれぞれx座標、y座標に読み替えたものの為、緯度経度位置の情報(分かりやすいのは標高とか)から作図しやすい反面、既存の地図に重ね合わせるには、図法の変換をせねばなりません。で、その変換をせねばならない用事が生まれましてね。やってみるわけです。

欲しいのは数式なんですけど、正距方位図法で検索してもなかなか該当する物が見つからず、azimuthal equidistant projectionで検索したら

Azimuthal Equidistant Projection [Wolfram MathWorld]

というドンピシャなページが引っかかりました。φとλはそれぞれ緯度、経度に相当する値で、それをx,yに変換したり、x,yから緯度・経度を求める式も載っています。すばらしいです。ということで、これを単純に実装するだけででき上がりました。

緯度・経度が-z座標、x座標になっているメッシュに対して適用します。
import math

app = Application
log = app.Logmessage

oObj = app.Selection(0)

pntPosArr = oObj.ActivePrimitive.Geometry.Points.PositionArray

newPosArr = [[],[],[]]

phi1 = XSIMath.DegreesToRadians(35.658056)
lamda0 =  XSIMath.DegreesToRadians(139.741389)

for i in range(len(pntPosArr[0])):
	phi = XSIMath.DegreesToRadians(-pntPosArr[2][i]) #latitude
	lamda = XSIMath.DegreesToRadians(pntPosArr[0][i]) #longitude
	
	c = math.acos( math.sin(phi1) * math.sin(phi) + math.cos(phi1) * math.cos(phi) * math.cos(lamda - lamda0) )
	
	k = c/math.sin(c)
	
	x = k * math.cos(phi) * math.sin(lamda - lamda0)
	z = k * (math.cos(phi1) * math.sin(phi) - math.sin(phi1) * math.cos(phi) * math.cos( lamda - lamda0 ))
	
	newPosArr[0].append( x )
	newPosArr[2].append( -z )
	newPosArr[1].append( 0.0 )
	
oObj.ActivePrimitive.Geometry.Points.PositionArray = newPosArr




上図の色がついている範囲の板(センターは原点にある)に適用すると



となる。見やすくするために明るさを持ち上げてあります。
形状も、冒頭にリンクを貼った日本周辺にほぼあっているかな、と思いますです、はい。
コメント ( 0 )|Trackback ( )
 
コメント
 
コメントはありません。
コメントを投稿する
ブログ作成者から承認されるまでコメントは反映されません
 
名前
タイトル
URL
コメント
コメント利用規約に同意の上コメント投稿を行ってください。

数字4桁を入力し、投稿ボタンを押してください。