# tessellate_plotΒΆ

This script plots the tessellated sphere into the 2-D map

To run this script, the Matplotlib Basemap module is needed to be installed.

conda install basemap    # Recommended


or

pip install https://github.com/matplotlib/basemap/archive/v1.1.0.tar.gz

#!/usr/bin/env python
''' This script plots the tessellated sphere into the 2-D map

To run this script, the Matplotlib Basemap module is needed to be installed.

::

conda install basemap    # Recommended

or

::

pip install https://github.com/matplotlib/basemap/archive/v1.1.0.tar.gz
'''

import sys
import time

import matplotlib
matplotlib.use('Agg')

import logging
logging.basicConfig()
logging.getLogger().setLevel(logging.DEBUG)

from pylab import *
import numpy as np

import matplotlib.patches as mpatches

from irfpy.util.tessellatedsphere import TessellatedSphere

def v3ang(v):
r = math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2])
t = math.acos(v[2] / r) * 180 / math.pi
p = math.atan2(v[1], v[0]) * 180 / math.pi

return (r, t, p)

def isnorm(x, y, z):
return (abs(x-y) <=180 and abs(y-z) < 180. and abs(z-x)< 180.)

def normalize_phi(p0, p1, p2):
''' Normalize if either or more p0=p1 p1=p2 p2=p0 is more than 180 deg.
'''

if isnorm(p0, p1, p2):
return p0, p1, p2
else:
p0m = p0 - 360
p0p = p0 + 360
p1m = p1 - 360
p1p = p1 + 360
p2m = p2 - 360
p2p = p2 + 360

if isnorm(p0m, p1, p2):
return (p0m, p1, p2)
elif isnorm(p0p, p1, p2):
return (p0p, p1, p2)
elif isnorm(p0, p1m, p2):
return (p0, p1m, p2)
elif isnorm(p0, p1p, p2):
return (p0, p1p, p2)
elif isnorm(p0, p1, p2m):
return (p0, p1, p2m)
elif isnorm(p0, p1, p2p):
return (p0, p1, p2p)

raise RuntimeError('Normalization failed.')

def plot_tesselation(ndim=0, filename=None):

logging.info('Given ndim=%d' % ndim)
tt0 = time.time()

tris = TessellatedSphere(ndim=ndim).getTriangles()

logging.info('Tessellation done.  %f sec' % (time.time() - tt0))

logging.info('Number of triangles is %d' % len(tris))

fig = figure()

from mpl_toolkits.basemap import Basemap
map = Basemap(projection='moll', lon_0=0.01, resolution = 'l', area_thresh=1000.)
map.drawmapboundary()

triarea = []

for tri in tris:
triarea.append(tri.getArea())

v0, v1, v2 = tri.getVertices()
r0, t0, p0 = v3ang(v0)   # r is radius, t is theta, p is phi
r1, t1, p1 = v3ang(v1)
r2, t2, p2 = v3ang(v2)

if abs(r0-1) > 1e-5 or abs(r1-1) > 1e-5 or abs(r2-1) > 1e-5:
raise RuntimeError('Failed! Not normalized! %s' % tri)

t0 = 90 - t0   # theta to be latitude.
t1 = 90 - t1
t2 = 90 - t2

lats = [t0, t1, t2, t0]
lons = [p0, p1, p2, p0]

if not isnorm(p0, p1, p2):
continue

try:
p0, p1, p2 = normalize_phi(p0, p1, p2)
except RuntimeError as e:
logging.warn('Failed the normalization in phi %s' % e)
continue

x, y = list(map(lons, lats))

map.plot(x, y, '-')

triarea = np.array(triarea)
print('Total area = %f' % sum(triarea))
print('Average area = %f' % average(triarea))
print('Total area of sphere = %f' % (4 * np.pi))
print('Average area of sphere with iso-area element = %f' % (4*np.pi / len(triarea)))

if filename != None:
fig.savefig(filename)
else:
plt.show()

if __name__ == '__main__':
if len(sys.argv) < 2:
print("USAGE:  %s   ndim   [figname.png]", file=sys.stderr)
sys.exit(-1)

if len(sys.argv) == 2:
filename = None
else:
filename = sys.argv[2]
plot_tesselation(ndim=int(sys.argv[1]), filename=filename)