tessellatedsphere_recplot

Recurrence version, plot.

''' Recurrence version, 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)

    ### 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))
    plt.savefig('tessellatedsphere_indecies_ndim{:02d}_{:d}x{:d}.png'.format(ndim, nlon, nlat))


def main():

    run2d(0)
    run2d(1)
    run2d(2)
    run2d(3)
    run2d(4)
    run2d(5)
    run2d(6)
    run2d(7)
    run2d(8)

    run2d(3, nlon=72, nlat=36)
    run2d(3, nlon=72*2, nlat=36*2)
    run2d(3, nlon=72*4, nlat=36*4)
    run2d(3, nlon=72*8, nlat=36*8)

    run(0)
    run(1)
    run(2)
#    run(3)


if __name__ == "__main__":
    main()