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)