apps130311_water.mapinme2

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 time
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, shuffle=False):
    logger = logging.getLogger(__name__)
    logger.setLevel(logging.DEBUG)

    ### Make directory
#    dirname = "mapinme_%dx%d_%s" % (nlon, nlat, datetime.datetime.now().strftime('%FT%T'))
    dirname = "mapinme2_%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()

    olist = list(range(orb0, orb1 + 1))
    if shuffle:
        import random
        random.shuffle(olist)

    for idone, orb in enumerate(olist):
        t1 = time.time()
        logger.info('Orbit %d/%d, Ellapsed=%f, Remains=%f' % (idone, len(olist), t1 - t0,
                        (t1 - t0) * len(olist) / (idone + 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))

            lonlist0 = (lonlist + dlon / 2.0)
            latlist0 = (latlist + dlat / 2.0)
            lon0d, lat0d = np.meshgrid(lonlist0, latlist0)

            lonlist0 = np.deg2rad(lonlist0)
            latlist0 = np.deg2rad(latlist0)
            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)

            matx_me2lse = irfpy.cy1orb.lseme.me2lse_matrix(tt)
            viewlse_full = matx_me2lse.dot(viewme.transpose((1, 0, 2)))   # (3, nlat, nlon)

            matx_lse2sc = irfpy.cy1orb.pvat2.get_lse2sc_matrix(tt)
            viewsc_full = matx_lse2sc.dot(viewlse_full.transpose((1, 0, 2)))   # (3, nlat, lon)

            matx_sc2cena = irfpy.cena.frame.get_sc2cena_matrix()
            viewcena_full = matx_sc2cena.dot(viewsc_full.transpose((1, 0, 2)))  # (3, nlat, nlon)

            theta_full, phi_full = irfpy.cena.frame.cena2angles2(
                            viewcena_full[0, :, :],
                            viewcena_full[1, :, :],
                            viewcena_full[2, :, :])
            theta_full = np.rad2deg(theta_full)   # (nlat, nlon)
            phi_full = np.rad2deg(phi_full)   # (nlat, nlon)

            theta_mask = ((theta_full <= -6.05 - 3.22) | (theta_full >= -6.05 + 3.22))
            phi_mask = ((phi_full <= -19.06 - 22.5) | (phi_full >= 19.06 + 22.5))
            angle_mask = angle.mask

            fullmask = theta_mask | phi_mask | angle_mask

            lonmask = np.ma.masked_array(lon0d, mask=fullmask).compressed()
            if lonmask.shape == (0,):
                continue
        
            latmask = np.ma.masked_array(lat0d, mask=fullmask).compressed()
            print(lonmask.shape, latmask.shape, file=sys.stderr)

#            for wlon, wlat in zip(lonmask, latmask):
#                m.append_value_lonlat(wlon, wlat, tt)
            ttarr = [tt for i in range(len(lonmask))]
            m.append_values_lonlat(lonmask, latmask, ttarr)

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



if __name__ == '__main__':
    t0 = time.time()
    r = 1
    runit(orb0=int(sys.argv[1]), orb1=int(sys.argv[2]), nlon=360 * r, nlat=180 * r, shuffle=True)
    print('Elapsed = %.2f' % (time.time() - t0))