apps130311_water.mapinme

Map all the orbit data to the surface.

Map all (cena mass mode data availble) the orbit data to the surface.

''' Map all the orbit data to the surface.

Map all (cena mass mode data availble) the orbit data to the surface.
'''

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
import irfpy.cy1orb.Cy1OrbitNr
from irfpy.cena import intersection

import irfpy.util.datafile
import irfpy.util.gridsphere

def masterinfo_reader(pathname):

    if os.path.isdir(pathname):
        pathname = os.path.join(pathname, 'masterinfo.dat')
    datf = irfpy.util.datafile.Datafile(open(pathname))
    datr = irfpy.util.datafile.DatafileReader(datf)
    return datr


def orbitdata_reader(pathname):
    fp = gzip.open(pathname)
    
    m = pickle.load(fp)   # SimpleGridSphere instance.

    return m


def runit(nlon=360, nlat=180, orb0=935, orb1=3160):
    logger = logging.getLogger(__name__)
    logger.setLevel(logging.DEBUG)

    ### Make directory
#    dirname = "mapinme_%dx%d_%s" % (nlon, nlat, datetime.datetime.now().strftime('%FT%T'))
    dirname = "mapinme_%dx%d" % (nlon, nlat)
    try:
        os.makedirs(dirname)
    except OSError as e:
        logger.warn('Some error, but continue anyway.\n\t(%s)' % e)
        pass

    # 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

    datf = irfpy.util.datafile.Datafile()
    datf.add_header('NLON', '%d' % nlon)
    datf.add_header('NLAT', '%d' % nlat)
    datf.add_header('ORBIT0', '%d' % orb0)
    datf.add_header('ORBIT1', '%d' % orb1)

    onr = irfpy.cy1orb.Cy1OrbitNr.Cy1OrbitNr()
    onr.setFromDefaultUri()
    starttime = onr.getStartTime(orb0)
    stoptime = onr.getStopTime(orb1)

    fp = open(os.path.join(dirname, 'masterinfo.dat'), 'w')
    datf.dump(fp)
    fp.close()

    t0 = time.time()

    for orb in range(orb0, orb1 + 1):
        t1 = time.time()
        logger.info('Orbit %d/%d, Ellapsed=%f, Remains=%f' % (orb - orb0, orb1 - orb0, t1 - t0,
                        (t1 - t0) * (orb1 - orb0) / (orb - orb0 + 1e-10)))
        m = irfpy.util.gridsphere.SimpleGridSphere(nlat=nlat, nlon=nlon)
    
        tlist = cm.getobstime(orbit=orb)

        if len(tlist) == 0:
            continue

        for tt in tlist:

            # 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))
            t2 = time.time()
            ttp = time.mktime(tt.timetuple())
            logger.debug('     Now at %s (%f). Ellpased %f s. Processed %f. Remains %f s' % (tt, ttp,
                                t2 - t0, (ttp - starttime) / (stoptime - starttime),
                                (t2 - t0) * (stoptime - ttp) / (ttp - starttime)))

            lonlist0 = (lonlist + dlon / 2.0) * np.pi / 180.
            latlist0 = (latlist + dlat / 2.0) * np.pi / 180.
            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)

            # 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

                    m.append_value_lonlat(np.rad2deg(lonrad), np.rad2deg(latrad), tt)

        fp = gzip.open(os.path.join(dirname, 'orbit%04d.out.gz' % orb), 'w')
        pickle.dump(m, fp)
        fp.close()



if __name__ == '__main__':
    import time
    t0 = time.time()
#    retval = main()
    runit(orb0=int(sys.argv[1]), orb1=int(sys.argv[2]))
    print('Elapsed = %.2f' % (time.time() - t0))