tess_plot_sample

''' A sample to plot tessellated data
'''
import time

import numpy as np
from matplotlib.mlab import griddata
import matplotlib.pyplot as plt
import matplotlib.cm

from irfpy.util import tessellatedsphere

def main():
    '''Main script'''

    ### Here you instance tessellated sphere of 4 dimension.
    ts = tessellatedsphere.RecurrentTessellatedUnitSphere(5)
    cog = ts.cogs    # (5120, 3) array
    loncog = np.arctan2(cog[:, 1], cog[:, 0]) * 180. / np.pi
    latcog = np.arcsin(cog[:, 2]) * 180. / np.pi

#    np.random.seed(100)
#    value = np.random.random(len(loncog))
#    print(value)

    _t0 = time.time()

    value = np.tanh(latcog / 45.) * np.cos(loncog / 12.)

    ### Grid data, the point of grid.

    dx = 2.5
    dy = 2.5
    xs, ys = np.mgrid[-90 + dx / 2.:90 + dx / 2. - dx / 1000.:dx,
                        0 + dy / 2.:50 + dy / 2. - dy / 1000.:dy]

    resampled = griddata(loncog, latcog, value, xs, ys)
#    print(resampled)

    plt.figure()
    plt.imshow(resampled.T, extent=(-90, 90, 0, 50), interpolation='nearest', origin='lower')
#    plt.scatter(loncog, latcog, edgecolor='none', color='r')
#    plt.scatter(xs, ys, color='b')
    plt.xlim(-92, 92)
    plt.ylim(-0.1, 50.1)

    _t1 = time.time()
    print('TIME I: {}'.format(_t1 - _t0))

    ### Another way... Manual sampling

    _t0 = time.time()

    # xs is lon, ys is lat, so the vector expression
    shape = xs.shape
    lonr = (xs * np.pi / 180.).flatten()
    latr = (ys * np.pi / 180.).flatten()
    vx = np.cos(lonr) * np.cos(latr)
    vy = np.sin(lonr) * np.cos(latr)
    vz = np.sin(latr)

    v = np.array([vx, vy, vz]).T    # (N, 3) array
    idx = ts.indecies(v)   # (N, array)
    value2 = value[idx]

    plt.figure()
    plt.imshow(value2.reshape(shape).T, extent=(-90, 90, 0, 50), interpolation='nearest', origin='lower')
#    plt.scatter(loncog, latcog, edgecolor='none', color='r')
#    plt.scatter(xs, ys, color='b')
    plt.xlim(-92, 92)
    plt.ylim(-0.1, 50.1)

    _t1 = time.time()
    print('TIME II: {}'.format(_t1 - _t0))

    ### Yet another way... Polygon

    from matplotlib.patches import Polygon
    from matplotlib.collections import PatchCollection

    _t0 = time.time()

    pcoll = []
    ccoll = []

    for idx, tri in enumerate(ts.tess_tris):
        vert = np.array(tri.getVertices())
        lonv = np.arctan2(vert[:, 1], vert[:, 0]) * 180. / np.pi
        latv = np.arcsin(vert[:, 2]) * 180. / np.pi

        lonlat = np.array([lonv, latv])    # (2, 3) array
#        print(lonlat)

        if lonv.max() - lonv.min() > 90:
            continue

        v = np.tanh(latv.mean() / 45.) * np.cos(lonv.mean() / 12.)
        ccoll.append(v)

        polygon = Polygon(lonlat.T, True)
        pcoll.append(polygon)

    p = PatchCollection(pcoll, cmap=matplotlib.cm.jet, edgecolor='none')
    p.set_array(np.array(ccoll))

    plt.figure()
    plt.gca().add_collection(p)

    plt.gca().set_aspect('equal')

    plt.xlim(-92, 92)
    plt.ylim(-0.1, 50.1)
#    plt.xlim(-180, 180)
#    plt.ylim(-90, 90)
    
    _t1 = time.time()
    print('TIME III: {}'.format(_t1 - _t0))


if __name__ == "__main__":
    main()