apps130311_water.pp

''' Simple plotting routine with polar projections
'''
import sys

import numpy as np
import matplotlib.pyplot as plt

import mapinmedata


def pol2xy_north(lonlist, latlist):
    ''' Projection. Polar (lonlat) to plot (XY).
    '''
    lonrarr = np.deg2rad(lonlist)
    colatarr = 90 - np.array(latlist)
    return colatarr * np.cos(lonrarr), colatarr * np.sin(lonrarr)

def test_pol2xy_north():
    lons = np.linspace(0, 360, 21)
    lats = np.linspace(-90, 90, 11)
    lonlist, latlist = np.meshgrid(lons, lats)
    lonlist = lonlist.flatten()
    latlist = latlist.flatten()
    x, y = pol2xy_north(lonlist, latlist)
    plt.figure()
    plt.scatter(x, y, c=lonlist)
    plt.colorbar()
    plt.title('LONGITUDE')
    plt.figure()
    plt.scatter(x, y, c=latlist)
    plt.colorbar()
    plt.title('LATITUDE')
        
def pol2xy_south(lonlist, latlist):
    ''' Projection. Polar (lo, lat) to plot (XY), in south pole.
    '''
    lonrarr = -np.deg2rad(lonlist)
    colatarr = 90 + np.array(latlist)
    return colatarr * np.cos(lonrarr), colatarr * np.sin(lonrarr)
        
def test_pol2xy_south():
    lons = np.linspace(0, 360, 21)
    lats = np.linspace(-90, 90, 11)
    lonlist, latlist = np.meshgrid(lons, lats)
    lonlist = lonlist.flatten()
    latlist = latlist.flatten()
    x, y = pol2xy_south(lonlist, latlist)
    plt.figure()
#    plt.scatter(x, y, c=lonlist)
    plt.contourf(x.reshape(11, 21), y.reshape(11, 21), lonlist.reshape(11, 21))
    plt.colorbar()
    plt.title('LONGITUDE')
    plt.figure()
#    plt.scatter(x, y, c=latlist)
    plt.contourf(x.reshape(11, 21), y.reshape(11, 21), latlist.reshape(11, 21))
    plt.colorbar()
    plt.title('LATITUDE')

def plot_northpole(lonarr, latarr, nlist, ax):

    xarr, yarr = pol2xy_north(lonarr, latarr)
    ax.contourf(xarr, yarr, np.ma.masked_less(nlist, 0.1))

    ## Latitude circle
    for lat in range(0, 90, 15):
        lon = np.linspace(0, 360, 90)
        xxx, yyy = pol2xy_north(lon, lat)
        ax.plot(xxx, yyy, 'k:')
    ## Longitude line
    for lon in range(0, 360, 30):
        lat = np.linspace(0, 90, 4)
        xxx, yyy = pol2xy_north(lon, lat)
        ax.plot(xxx, yyy, 'k:')

    ## Boundary lat = 60 deg here.
    ax.set_xlim(-60, 60)
    ax.set_ylim(-60, 60)
    ax.set_aspect(1)
    return

def plot_southpole(lonarr, latarr, nlist, ax):
    xarr, yarr = pol2xy_south(lonarr, latarr)
    ax.contourf(xarr, yarr, np.ma.masked_less(nlist, 0.1))

    ## Latitude line
    for lat in range(0, 90, 15):
        lon = np.linspace(0, 360, 90)
        xxx, yyy = pol2xy_north(lon, lat)
        ax.plot(xxx, yyy, 'k:')
    ## Longitude line
    for lon in range(0, 360, 30):
        lat = np.linspace(0, 90, 4)
        xxx, yyy = pol2xy_north(lon, lat)
        ax.plot(xxx, yyy, 'k:')

    ax.set_xlim(-60, 60)
    ax.set_ylim(-60, 60)
    ax.set_aspect(1)
    return

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

    m = mapinmedata.orbitdata_reader(sys.argv[1])
    lonarr, latarr = m.get_cgrid()
    nlist = m.get_counts()

    fig = plt.figure(figsize=(12, 6))
    ax1 = fig.add_subplot(1, 2, 1)
    ax2 = fig.add_subplot(1, 2, 2)

    plot_northpole(lonarr, latarr, nlist, ax1)
    plot_southpole(lonarr, latarr, nlist, ax2)
    plt.savefig(sys.argv[1]+'_pc.png')


if __name__ == "__main__":
    main()