from mpl_toolkits.basemap import Basemap, _geoslib
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
import numpy as np
from numpy import pi
icosahedron = \
[[0.53,0.,-0.53,0.53,-0.53,0.,0.53,-0.53,0.,0.53,0.,
-0.53,0.85,0.53,0.85,0.85,0.85,0.53,-0.85,-0.53,-0.85,
-0.85,-0.85,-0.53,0.,0.85,0.,0.,0.,-0.85,0.,0.,0.85,
0.,-0.85,0.,0.,0.53,0.85,0.,0.85,0.53,0.53,0.,0.85,
0.53,0.85,0.,-0.53,0.,-0.85,-0.53,-0.85,0.,-0.53,
-0.85,0.,-0.53,0.,-0.85],
[0.,0.85,0.,0.,0.,-0.85,0.,0.,0.85,0.,-0.85,0.,0.53,
0.,-0.53,0.53,-0.53,0.,-0.53,0.,0.53,-0.53,0.53,0.,
0.85,0.53,0.85,0.85,0.85,0.53,-0.85,-0.85,-0.53,-0.85,
-0.53,-0.85,0.85,0.,0.53,0.85,0.53,0.,0.,-0.85,-0.53,
0.,-0.53,-0.85,0.,0.85,0.53,0.,0.53,0.85,0.,-0.53,
-0.85,0.,-0.85,-0.53],
[0.85,0.53,0.85,0.85,0.85,0.53,-0.85,-0.85,-0.53,-0.85,
-0.53,-0.85,0.,0.85,0.,0.,0.,-0.85,0.,0.85,0.,0.,0.,
-0.85,0.53,0.,-0.53,0.53,-0.53,0.,0.53,-0.53,0.,0.53,
0.,-0.53,0.53,0.85,0.,-0.53,0.,-0.85,0.85,0.53,0.,
-0.85,0.,-0.53,0.85,0.53,0.,-0.85,0.,-0.53,0.85,0.,
0.53,-0.85,-0.53,0.]]
icosahedron1 = \
[[0.53, 0. ,-0.53, 0.53,-0.53, 0. ],
[0. , 0.85, 0. ,0. , 0. ,-0.85],
[0.85, 0.53, 0.85 ,0.85, 0.85, 0.53]]
def pointsOnSphere():
x,y,z = np.array(icosahedron)/0.851
return 180/pi*np.arcsin(z), 180/pi*np.arctan2(y,x)
if __name__=='__main__':
if 0:
from enthought.mayavi import mlab
x,y,z = icosahedron
sphere = mlab.triangular_mesh(x, y, z, \
np.arange(len(x)).reshape(-1,3), representation = 'wireframe')
mlab.show()
raise SystemExit
# Make Mollweide plot
m = Basemap(projection='moll', lon_0=0, resolution='c')
# draw the edge of the map projection region (the projection limb)
m.drawmapboundary()
theta, phi = pointsOnSphere()
theta.shape = phi.shape = (-1,3)
print theta.shape[0], 'polys'
ax = plt.gca() # get current axes instance
# create a geos Polygon instance describing the map region
boundarypoly = _geoslib.Polygon(m._boundarypolyxy.boundary)
for i in range(theta.shape[0]):
pts = np.vstack((phi[i], theta[i])).T
polypts = np.array([m(*pt) for pt in pts]) # to projection coords
poly = _geoslib.Polygon(polypts) # geos polygon for testing
if poly.intersects(boundarypoly):
for psub in poly.intersection(boundarypoly):
b = psub.boundary # boundary of an intersection
polypts = zip(b[:,0],b[:,1])
c = (1,) + tuple(np.random.random(2)) # warm colour
poly = Polygon(polypts, facecolor=c, edgecolor=c)
ax.add_patch(poly)
else:
c = tuple(np.random.random(2)) + (1,) # cool colour
poly = Polygon(polypts, facecolor=c, edgecolor=c)
ax.add_patch(poly)
plt.show()