python で球面グリッド三角化

```# coding: utf-8

# In[333]:

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import pyshtools as sh
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import Delaunay
import math
import collections # 多値返却のため
get_ipython().magic(u'matplotlib notebook')

# In[334]:

# 経線、緯線数を与えて三角グリッドを作る
def spheregrid(nlat=10,nlon=10):
phi = 2*math.pi/nlon * np.array([range(nlon)])
theta = math.pi/(nlat-1) * np.array([range(nlat)])

n = 2 + nlon * (nlat-2)

ret = np.zeros((2,n))
ret[0,n-1] = math.pi
for i in range(nlat -2):
for j in range(nlon):
tmp = i * nlon + j + 1
ret[0,tmp] = theta[0,i+1]
ret[1,tmp] = phi[0,j]
n_tri = nlon * 2 + nlon*2*(nlat-3)
tri = np.zeros((n_tri,3),dtype=np.int)
cnt = 0
for i in range(nlon-1):
tri[cnt,0] = 0
tri[cnt,1] = i+1
tri[cnt,2] = i+2
cnt +=1
tri[cnt,0] = 0
tri[cnt,1] = nlon
tri[cnt,2] = 1
cnt +=1
for i in range(nlat-3):
for j in range(nlon-1):
tri[cnt,0] = i*nlon + j + 1
tri[cnt,1] = i*nlon + j + 2
tri[cnt,2] = (i+1)*nlon + j + 2
cnt += 1
tri[cnt,0] = i*nlon + j + 1
tri[cnt,1] = (i+1)*nlon + j + 1
tri[cnt,2] = (i+1)*nlon + j + 2
cnt += 1
tri[cnt,0] = i*nlon + nlon
tri[cnt,1] = i*nlon + 1
tri[cnt,2] = (i+1)*nlon + 1
cnt += 1
tri[cnt,0] = i*nlon + nlon
tri[cnt,1] = (i+1)*nlon + nlon
tri[cnt,2] = (i+1)*nlon + 1
cnt += 1
for i in range(nlon-1):
tri[cnt,0] = (nlat-3)*nlon + i + 1
tri[cnt,1] = (nlat-3)*nlon + i + 2
tri[cnt,2] = n-1
cnt +=1
tri[cnt,0] = (nlat-3)*nlon + nlon-1 + 1
tri[cnt,1] = (nlat-3)*nlon + 1
tri[cnt,2] = n-1

X = np.ndarray((3,n))
for i in range(ret[0,].size):
X[0,i], X[1,i], X[2,i] = np.sin(ret[0,i])*np.cos(ret[1,i]), np.sin(ret[0,i])*np.sin(ret[1,i]), np.cos(ret[0,i])

result = collections.namedtuple('result', 'X, angles, tri')
return result(X=X, angles=ret, tri=tri)

# In[335]:

# 3d 正規乱点を単位球面に発生させる
def runitsphere(n,d=3):
X = np.random.randn(n*d).reshape(n,d)
r = np.sum(np.power(X,2), axis=1)
X = X.T/np.power(r,0.5)
return X

# In[336]:

# 球面上乱点を発生しその3角化をする関数を作る
def spheretri_random(n,s=0.001):
X = runitsphere(n-1) # 点 (0,0,1)を加えて全部でn個の点を球面に置くためにn-1とする
# (0,0,1)のごく近辺に3点を発生させる
X_min = np.min(np.abs(X))
xy_min = X_min * s
three1x = xy_min
three1y = 0
three1z = math.sqrt(1-three1x**2-three1y**2)
three2x = xy_min * (-1/2)
three2y = xy_min * (math.sqrt(3)/2)
three2z = math.sqrt(1-three2x**2-three2y**2)
three3x = xy_min * (-1/2)
three3y = xy_min * (-math.sqrt(3)/2)
three3z = math.sqrt(1-three3x**2-three3y**2)
# 2通りの球面上点座標を作る
# (0,0,1)に近い3点をそれぞれ登録する場合
x2 = np.append(X[0,],np.array([three1x,three2x,three3x]))
y2 = np.append(X[1,],np.array([three1y,three2y,three3y]))
z2 = np.append(X[2,],np.array([three1z,three2z,three3z]))
# 3点を登録する代わりに(0,0,1)を登録する場合
x = np.append(X[0,],np.array([0]))
y = np.append(X[1,],np.array([0]))
z = np.append(X[2,],np.array([1]))
# (0,0,1)近傍3点を区別する場合に関して
# (0,0,1) が無限遠となるように、球面を平面に開く
# そうすることでその二次元平面にて三角化する
# その上で(0,0,1)相当の3点は、最終的に、三角を作るものとする
X2 = 2*x2/(2-(z2+1))
Y2 = 2*y2/(2-(z2+1))
# ドロネー3角化(二次元平面にて行う)
# Triangulate parameter space to determine the triangles
tri = Delaunay(np.array([X2,Y2]).T)
# 最後に追加した3点は、(0,0,1)のことなので、3角化の頂点座標のうち、追加点のノードIDを同じ値にする
tri.simplices[tri.simplices > n-2]=n-1
# 返り値は、n個の点の座標とその三角形頂点リスト
result = collections.namedtuple('result', 'X, tri')
retX = np.vstack([np.vstack([x,y]),z])
return result(X=retX, tri=tri.simplices)

# 三角形の面積を計算する
def triarea(v1,v2,v3):
e12 = v2-v1
e13 = v3-v1
ip = np.dot(e12,e13)
r1 = np.sqrt(np.sum(np.power(e12,2)))
r2 = np.sqrt(np.sum(np.power(e13,2)))
costheta = ip/(r1*r2)
sintheta = np.sqrt(1-np.power(costheta,2))
area = 0.5 * r1*r2*sintheta
return area

# 三角形の面積を三角化オブジェクトから作る
def triarea_ts(ts):
ret = np.zeros(ts.tri.T[0,].size)
for i in range(ts.tri.T[0,].size):
id1 = ts.tri[i,0]
id2 = ts.tri[i,1]
id3 = ts.tri[i,2]
if (id1-id2)*(id2-id3)*(id3-id1) == 0:
ret[i] = 0
else:
v1 = np.array([ts.X[0,id1],ts.X[1,id1],ts.X[2,id1]])
v2 = np.array([ts.X[0,id2],ts.X[1,id2],ts.X[2,id2]])
v3 = np.array([ts.X[0,id3],ts.X[1,id3],ts.X[2,id3]])
ret[i] = triarea(v1,v2,v3)

return ret

# 三角形の重心座標を三角化オブジェクトから作る
def triweightedcenter(ts):
ret = np.zeros((3,ts.tri.T[0,].size))
for i in range(ts.tri.T[0,].size):
id1 = ts.tri[i,0]
id2 = ts.tri[i,1]
id3 = ts.tri[i,2]

v1 = np.array([ts.X[0,id1],ts.X[1,id1],ts.X[2,id1]])
v2 = np.array([ts.X[0,id2],ts.X[1,id2],ts.X[2,id2]])
v3 = np.array([ts.X[0,id3],ts.X[1,id3],ts.X[2,id3]])
ret[0,i], ret[1,i], ret[2,i] = (v1 + v2 + v3)/3

return ret

# 三角化プロット
def spheretriplot(ts):
# 3Dplot
fig = plt.figure()
ax = Axes3D(fig)
# The triangles in parameter space determine which x, y, z points are
# connected by an edge
ax.plot_trisurf(ts.X[0,],ts.X[1,],ts.X[2,], triangles=ts.tri, cmap=plt.cm.Spectral)
#ax.plot_trisurf(ts.X[0,],ts.X[1,],ts.X[2,], triangles=ts.tri, color=col)
#ax.plot_trisurf(x, y, z, triangles=tri.simplices)
plt.show()

# In[337]:

# 3Dデカルト座標を極座標に変える(r==0は考えない)
def decart2polar3d(x,y,z):
r = math.sqrt(x**2+y**2+z**2)
rxy = x**2+y**2
if rxy==0:
phi = 0
theta = 0
if z < 0:
theta = math.pi
else:
if x==0:
phi = math.pi/2
if y < 0:
phi = -math.pi/2
else:
phi = np.arctan(y/x)
theta = np.arccos(z/r)
return [phi,theta,r]

decart2polar3d(0,0,-1)

def decart2polar3d_ndarray(X):
ret = np.zeros_like(X)
for i in range(X[0,].size):
tmp= decart2polar3d(X[0,i],X[1,i],X[2,i])
ret[0,i],ret[1,i],ret[2,i] = tmp
return ret

# In[338]:

# 三角の頂点・重心座標の場の値を球面調和関数係数numpy.ndarrayから作る
def SHval(Xpolar,coef):
valreal = np.zeros(Xpolar[0,].size)
valcomplex = np.zeros(Xpolar[0,].size,'complex')
for i in range(Xpolar[0,].size):
#valreal += [ coeff2valreal(coef,Xpolar[0,i],Xpolar[1,i])]
#valcomplex += [coeff2valcomplex(coef,Xpolar[0,i],Xpolar[1,i])]
valreal[i] = coeff2valreal(coef,Xpolar[0,i],Xpolar[1,i])
valcomplex[i] = coeff2valcomplex(coef,Xpolar[0,i],Xpolar[1,i])

result = collections.namedtuple('result', 'real, clx')
return result(real = valreal,clx=valcomplex)

# In[339]:

# 球面になめらかな場を与える
# 球面調和関数の係数のセットとして与える
# x.coeffsは[0,,]にm>=0[1,,]にm<0が格納される
# [*,0,k]はabs(m)=kの値を格納する
def randSHcoef(lmax):
num_coef = (lmax+1)**2
value_re = np.random.randn(num_coef)
#print value_re
x = sh.SHCoeffs.from_zeros(lmax)

for l in range(0, lmax+1):
for m in range(-l, l+1):
x.set_coeffs(value_re[l*(l+1)+m], l, m)
return x

# In[340]:

# scipy.special.sph_harmを使ってSHcoeffから、球面座標の値を計算
def coeff2valreal(cf,phi,theta):
ret = 0
#lmax = cf
for l in range(0, cf[0,0,].size):
for m in range(-l, l+1):
if m >= 0:
tmp = cf[0,l,m]
ret += sph_harm(m,l,phi,theta).real
else:
tmp = cf[1,l,np.abs(m)]
ret += sph_harm(m,l,phi,theta).imag
return ret

def coeff2valcomplex(cf,phi,theta):
ret = 0+1j
#lmax = cf
for l in range(0, cf[0,0,].size):
for m in range(-l, l+1):
if m >= 0:
tmp = cf[0,l,m]
else:
tmp = cf[1,l,np.abs(m)]
ret += sph_harm(m,l,phi,theta)
return ret

# # 球面上一様乱点とその三角化オブジェクト tsを作る
# n = 200
# ts2 = spheretri_random(n)
# #print ts2
# ts = spheregrid(nlat=10,nlon=10)
# #print ts
# #print(ts.tri)
# #print ts.tri.T[0,].size
# # 三角形の面積
# areas = triarea_ts(ts)
# # 三角形の重心座標
# weightedcenter = triweightedcenter(ts)
# #print weightedcenter
# # 頂点と三角形重心の極座標
# Xpolar = decart2polar3d_ndarray(ts.X)
# #Xpolar = ts.angles
# WCpolar = decart2polar3d_ndarray(weightedcenter)
# #print weightedcenter.size
# #print WCpolar.size
# spheretriplot(ts)

# In[343]:

# 球面場モデルを作る
f = randSHcoef(5)

# In[344]:

# 場の値を計算
fvalV = SHval(Xpolar,f.coeffs)
fvalF = SHval(WCpolar,f.coeffs)
#print WCpolar[0,].size
#print WCpolar
#print fvalV
#print fvalF.clx.size

# In[ ]:

# In[345]:

x = randSHcoef(3)
print x.coeffs
# 係数はそのまま回転できる
x_rot = x.rotate(1,2,3, degrees=False)
print x_rot.coeffs

xexp = x.expand()

# In[346]:

coef = randSHcoef(3).coeffs
#coef[0,0,].size
print coeff2valreal(coef,Xpolar[0,0],Xpolar[1,0])
print coeff2valcomplex(coef,Xpolar[0,0],Xpolar[1,0])
#Xpolar[0,0]

valreal = np.zeros(Xpolar[0,].size)
valcomplex = np.zeros(Xpolar[0,].size,'complex')
for i in range(Xpolar[0,].size):
#valreal += [ coeff2valreal(coef,Xpolar[0,i],Xpolar[1,i])]
#valcomplex += [coeff2valcomplex(coef,Xpolar[0,i],Xpolar[1,i])]
valreal[i] = coeff2valreal(coef,Xpolar[0,i],Xpolar[1,i])
valcomplex[i] = coeff2valcomplex(coef,Xpolar[0,i],Xpolar[1,i])

# In[348]:

plt.figure()
plt.plot(valreal,valcomplex.real)
plt.show()
plt.figure()
plt.plot(valcomplex.real,valcomplex.imag)
plt.show()

# In[ ]:

```