tessellatedsphere_pyplot

""" Example of plotting.
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
from irfpy.util import nptools
import matplotlib.pyplot as plt
import os, sys
os.environ['PROJ_LIB'] = os.path.join(sys.exec_prefix, 'share', 'proj')
from mpl_toolkits.basemap import Basemap
from irfpy.util.tessellatedsphere import RecurrentTessellatedUnitSphere as RTsphere

def main():

    ts = RTsphere(ndim=6)
    lonlist, latlist = ts.getTriangleLonLat()
    n = lonlist.shape[0]
    print('n={}'.format(n))
    #vals = np.arange(n)    # Value is the triangle ID
    vals = np.random.uniform(size=n)    # Value is the triangle ID
    
    # plot vertex on usual matplotlib
    
    plt.scatter(lonlist, latlist, s=1)
    plt.savefig('tessellatedsphere_pyplot_fig1.png')
    plt.close('all')
    
    # Using triplot on usual matplotlib
    
    x_coords_full = lonlist.flatten()
    y_coords_full = latlist.flatten()
    tris_full = np.arange(n * 3).reshape(n, 3)
    
    ## Here some tris shall be removed...
    x_coords = []
    y_coords = []
    tris_valid = []
    vals_valid = []
    
    for i in range(n):
        _lon = nptools.unjump_periodic(lonlist[i], 180)   # (3,)
        _lat = latlist[i]   # (3,)
        _tris = len(x_coords)
    #    _vals = i
        _vals = vals[i]
    
        if _lon.max() - _lon.min() <= 180:   # No hemispheric transverse
            x_coords.append(_lon[0])
            x_coords.append(_lon[1])
            x_coords.append(_lon[2])
            y_coords.append(_lat[0])
            y_coords.append(_lat[1])
            y_coords.append(_lat[2])
            tris_valid.append([_tris, _tris+1, _tris+2])
            vals_valid.append(_vals)
        else:
            pass
            
    
    plt.tripcolor(np.array(x_coords), np.array(y_coords), np.array(tris_valid), np.array(vals_valid))
    # Trick here is to plot this two more times to consider the periodicity...
    plt.tripcolor(np.array(x_coords) - 360, np.array(y_coords), np.array(tris_valid), np.array(vals_valid))
    plt.tripcolor(np.array(x_coords) + 360, np.array(y_coords), np.array(tris_valid), np.array(vals_valid))
    plt.xlim(-180, 180)
    plt.ylim(-90, 90)
    plt.savefig('tessellatedsphere_pyplot_fig2.png')
    plt.close('all')
    
    # Global Basemap, vertex
    

    m = Basemap(projection='hammer', lon_0=0, resolution='c')
    m.drawcoastlines()
    
    mlon, mlat = m(lonlist, latlist)
    
    m.scatter(mlon, mlat, s=1)
    
    plt.savefig('tessellatedsphere_pyplot_fig3.png')
    plt.close('all')
    
    # Global Basemap, tripcolor
    m = Basemap(projection='nplaea', boundinglat=70, lon_0=0, resolution='c')
    m.drawcoastlines()
    
    mlon, mlat = m(lonlist, latlist)
    print(mlon.shape)
    
    x_coords = []
    y_coords = []
    tris_valid = []
    vals_valid = []
    for i in range(mlon.shape[0]):
    #    print(mlon[i], mlat[i])
        if np.isfinite(mlon[i]).all() and np.isfinite(mlat[i]).all():
            _lon = mlon[i]
            _lat = mlat[i]
            _tris = len(x_coords)
            _vals = vals[i]
    
    #        if _lon.max() - _lon.min() <= 180:   # No hemispheric transverse
            x_coords.append(_lon[0])
            x_coords.append(_lon[1])
            x_coords.append(_lon[2])
            y_coords.append(_lat[0])
            y_coords.append(_lat[1])
            y_coords.append(_lat[2])
            tris_valid.append([_tris, _tris+1, _tris+2])
            vals_valid.append(_vals)
        else:
            print('SKIP', i, mlon[i], mlat[i], )
    
    plt.tripcolor(np.array(x_coords), np.array(y_coords), np.array(tris_valid), np.array(vals_valid))
    plt.colorbar()
    
    plt.savefig('tessellatedsphere_pyplot_fig4.png')
    
    plt.close('all')
    
    
if __name__ == "__main__":
    main()