apps130311_water.esoa_1_3_number_finitefov_in_grid

Examine number of projection.

Output is a graphics ‘esoa_1_3_number_finitefov_in_grid.png’ and data file ‘esoa_1_3_number_finitefov_in_grid.out.gz’.

Output data file can be loaded via esoalib.esoa_1_3_load function.

''' Examine number of projection.

Output is a graphics 'esoa_1_3_number_finitefov_in_grid.png' and
data file 'esoa_1_3_number_finitefov_in_grid.out.gz'.

Output data file can be loaded via esoalib.esoa_1_3_load function.
'''

module_debug = True

import os
import sys
import logging
logging.basicConfig()
import datetime
import math

import gzip
import pickle

import bisect

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

import irfpy.cena.cena_mass2 as cm
import irfpy.cy1orb.pvat2 as pvat
import irfpy.cy1orb.lseme
from irfpy.cena import intersection


def runit(orbs, nlon=360, nlat=180):
    logger = logging.getLogger(__name__)
    logger.setLevel(logging.DEBUG)

    # Cartesian grid preparation

    # each list has the "lowest" value of each grid.
    lonlist = np.linspace(0, 360, num=nlon, endpoint=False)
    latlist = np.linspace(-90, 90, num=nlat, endpoint=False)
    dlon = 360.0 / nlon
    dlat = 180.0 / nlat

    logger.debug('LON=%s, (dLon=%.2f)' % (lonlist, dlon))
    logger.debug('LAT=%s, (dLat=%.2f)' % (latlist, dlat))

    ndata_grid = np.zeros([nlon, nlat])

    lonlist0 = (lonlist + dlon / 2.0) * np.pi / 180.
    latlist0 = (latlist + dlat / 2.0) * np.pi / 180.

    # Prepare a list
    time_position_data = [[[] for ilat in range(nlat)] for ilon in range(nlon)]

    for orb in orbs:
        tlist = cm.getobstime(orbit=orb)
        for tt in tlist[500:530:3]:  # Only 10 scans.
#        for tt in tlist:
            if module_debug:
                import warnings
                warnings.warn('MODULE DEBUG')

            print(tt)
            
#        for tt in tlist:
            # Only for SZA > 80 is taken
            sza = pvat.getsza(tt)
            if sza > 80:
                continue

            # Position of the S/C in ME frame as a vector
            sclon, sclat, schei = pvat.getmepos(tt)
            scposme = pvat.getmepos(tt, as_vector=True)
            # height is
            height = np.sqrt((scposme ** 2).sum())

            # Horizon angle is calculated.
            horizon_angle = np.arccos(1738. / height)
            logger.debug('Lon=%f, Lat=%f, Height = %f, Horizon Angle = %f' % (sclon, sclat, height, horizon_angle * 180. / np.pi))

            lon0, lat0 = np.meshgrid(lonlist0, latlist0)  # (nlat, nlon) shape.

            # vsurf is the vector from the center of the Moon to a grid
            vsurf = 1738. * np.array([np.cos(lat0) * np.cos(lon0), np.cos(lat0) * np.sin(lon0), np.sin(lat0)])    # shape of (3, nlat, nlon)
            vsurflen = np.sqrt((vsurf * vsurf).sum(0))   # (nlat, nlon) shape


            angle = np.arccos(vsurf.T.dot(scposme).T / height / vsurflen)   # (nlat, nlon)
            angle = np.ma.masked_greater(angle, horizon_angle)    # (nlat, nlon)


            viewme = (vsurf.T - scposme).T   # (3, nlat, nlon)

            # Confirmed by here.
#            from pudb import set_trace; set_trace()

            # Loop in (LON LAT) spece
            for ilon, lonrad in enumerate(lonlist0):
                for ilat, latrad in enumerate(latlist0):

                    if angle.mask[ilat, ilon]:
                        continue

                    # ME to lse
                    viewlse = irfpy.cy1orb.lseme.me2lse(tt, viewme[:, ilat, ilon])
                    # LSE to SC
                    viewsc = irfpy.cy1orb.pvat2.getscvec(tt, viewlse)
                    # SC to cena
                    viewcena = irfpy.cena.frame.sc2cena(viewsc)
                    theta, phi = irfpy.cena.frame.cena2angles(viewcena, degrees=True)

                    if theta <= -6.05 - 3.22 or theta >= -6.05 + 3.22:
                        continue
                    if phi <= -19.06 - 22.5 or phi >= 19.06 + 22.5:
                        continue

                    ndata_grid[ilon, ilat] += 1

                    # Save the time to the corresponding longitude and latitude
                    time_position_data[ilon][ilat].append(tt)

#
#            # A first step: spacecraft position
#            # lon is -180 to 180, lat is -90 to 90.
#            lon, lat = intersection.get_center_me(tt, 3)
#            while lon >= 360:
#                lon -= 360
#            while lon < 0:
#                lon += 360
#
#            # Then, map to grid
#            ilon = bisect.bisect(lonlist, lon) - 1
#            ilat = bisect.bisect(latlist, lat) - 1
#
#            logger.debug([tt, lon, lat, ilon, ilat, lonlist[ilon], latlist[ilat]])
#            
#            ndata_grid[ilon, ilat] += 1

    return lonlist, latlist, ndata_grid, time_position_data


def main():

    orbs = [2646,]

#    orbs = [2646,]
    nlon = 360 * 1
    nlat = 180 * 1
#    nlon = 90
#    nlat = 45

    lonlist, latlist, ndata_grid, time_position_data = runit(orbs, nlon=nlon, nlat=nlat)

    lonlist = list(lonlist)
    lonlist.append(360)

    latlist = list(latlist)
    latlist.append(90)

    x, y = np.meshgrid(lonlist, latlist)

    ndata_grid = np.ma.masked_equal(ndata_grid, 0)

    fig = plt.figure()
    ax = fig.add_subplot(111)
    pc = ax.pcolor(x, y, ndata_grid.T)

    ax.set_xlim(0, 360)
#    ax.set_xlim(220, 280)
    ax.set_ylim(-90, 90)
#    ax.set_ylim(-60, -30)
    ax.grid()

    clb = fig.colorbar(pc)
    clb.set_label('# CENA data (mapped for FOV CH-234)')

    fig.savefig('./esoa_1_3_number_finitefov_in_grid.png')

    fp = gzip.open('./esoa_1_3_number_finitefov_in_grid.out.gz', 'w')

    pickle.dump(time_position_data, fp)
    pickle.dump(ndata_grid, fp)
    pickle.dump(x, fp)  # lonlist for pcolor
    pickle.dump(y, fp)  # latlist for pcolor
    fp.close()

    return lonlist, latlist, ndata_grid, time_position_data


if __name__ == '__main__':
    import time
    t0 = time.time()
    retval = main()
    print('Elapsed = %.2f' % (time.time() - t0))