tessellatedsphere_conplot

Recurrence version, connection plot.

''' Recurrence version, connection plot.
'''
import time
import matplotlib.pyplot as plt
import numpy as np

from irfpy.util import tessellatedsphere

def run(ndim):
    tn = tessellatedsphere.RecurrentTessellatedUnitSphere(ndim=ndim)

    ### Grid.
    longrid = np.linspace(0, 360, 72)
    latgrid = np.linspace(-90, 90, 36)

    longridc = (longrid[1:] + longrid[:-1]) / 2.
    latgridc = (latgrid[1:] + latgrid[:-1]) / 2.

    lon, lat = np.meshgrid(longridc, latgridc)
    lonr, latr = np.deg2rad((lon, lat))

    gridshape = lon.shape
    ngrid = gridshape[0] * gridshape[1]

    xs = np.cos(lonr) * np.cos(latr)
    ys = np.sin(lonr) * np.cos(latr)
    zs = np.sin(latr)

    idxn = np.zeros(gridshape)
    
    _t0 = time.time()
    for i in range(gridshape[0]):
      for j in range(gridshape[1]):
        idxn[i, j] = tn.index1(np.array([xs[i, j], ys[i, j], zs[i, j]]))
    _t1 = time.time() - _t0
    print('time (ndim={}; index1) = {:f}s.  (1 search={:f}s)'.format(ndim, _t1, _t1 / ngrid))

    plt.pcolor(longrid, latgrid, idxn)
    plt.savefig('tessellatedsphere_index1_ndim{:02d}.png'.format(ndim))

def run2d(ndim, nlon=72, nlat=36):
    tn = tessellatedsphere.RecurrentTessellatedUnitSphere(ndim=ndim)
    tn.createLinkList()

    ### Grid.
    longrid = np.linspace(0, 360, nlon)
    latgrid = np.linspace(-90, 90, nlat)

    longridc = (longrid[1:] + longrid[:-1]) / 2.
    latgridc = (latgrid[1:] + latgrid[:-1]) / 2.

    lon, lat = np.meshgrid(longridc, latgridc)
    lonr, latr = np.deg2rad((lon, lat))

    gridshape = lon.shape
    ngrid = gridshape[0] * gridshape[1]

    lonr = lonr.flatten()
    latr = latr.flatten()

    xs = np.cos(lonr) * np.cos(latr)
    ys = np.sin(lonr) * np.cos(latr)
    zs = np.sin(latr)

    points = np.array([xs, ys, zs]).T   # (ngrid, 3) array
#    print(points.shape)
    
    _t0 = time.time()
    idxn = tn.indecies(points)
    _t1 = time.time() - _t0
    print('time (ndim={}; indecies, {nlon}x{nlat}) = {:f}s.  (1 search={:f}s)'.format(ndim, _t1, _t1 / ngrid, nlon=nlon, nlat=nlat))

    plt.pcolor(longrid, latgrid, idxn.reshape(gridshape))
    conn = tn.linklist    # (ntri, 3) array

    print(tn.cogs.shape)

    ntri = conn.shape[0]
    for itri in range(ntri):
        gc = tn.cogs[itri]
        lon = np.arctan2(gc[1], gc[0]) * 180. / np.pi
        lat = np.arcsin(gc[2] / np.sqrt((gc ** 2).sum())) * 180. / np.pi
        while lon < 0:
            lon += 360

        for dest in range(3):
            gc = tn.cogs[conn[itri, dest]]
            lond = np.arctan2(gc[1], gc[0]) * 180. / np.pi
            latd = np.arcsin(gc[2] / np.sqrt((gc ** 2).sum())) * 180. / np.pi
            while lond < 0:
                lond += 360

            print(lon, lond)
            print(lat, latd)
            if np.abs(lon - lond) <= 180:
                plt.plot([lon, lond], [lat, latd], '-')


    plt.savefig('tessellatedsphere_con_ndim{:02d}.png'.format(ndim))


def main():

    run2d(1)
    run2d(2)


if __name__ == "__main__":
    main()