'''Chandrayaan-1 orbital position
This will provide low-level access of orbit/attitude data.
For practical use, you may use high-level module, :mod:`irfpy.cy1orb.pvat2`.
.. note::
This module should *not* be used for data analysis.
Use :mod:`irfpy.cy1orb.pvat2`
See also :ref:`get_start_positon`.
The :class:`Cy1PvatRaw` and :class:`Cy1Pvat` classes are the lowest-API classes.
They represent a single master orbit-attitude file.
The master files are generated by some other tools (sorry, I forgot the software...
We have dataset!).
:class:`Cy1Orbit` and :class:`Cy1Attitude` is a second lowest-API classes.
'''
import os
import sys
import logging
import urllib.request
import urllib.parse
import urllib.error
import numpy
import dateutil.parser
import gzip
import pickle as pickle
from irfpy.util.irfpyrc import Rc
from irfpy.util.julday import JdSeries, JdObject, Julday
from irfpy.util.vector3d import Vector3d
[docs]class Cy1PvatRaw:
'''CY-1 PVAT file reader.
PVAT file should be a format of 23-column ascii file that SARA spice module will produce.
The data is normally read through :meth:`getPvatOfOrbit` method:
>>> pvat=Cy1PvatRaw()
>>> data=pvat.getPvatOfOrbit(1234)
In this case, orbit 1234 data is loaded into Ntime x 22 matrix.
One column (the 3rd column) is disregarded because it is not float.
In the above case, data is stored into the np.array instance.
The array is cached into memory in the instance of :class:`Cy1PvatRaw`.
:meth:`clearCache` method will remove the memory cache.
>>> pvat = Cy1PvatRaw()
>>> onr = 1280
>>> data = pvat.getPvatOfOrbit(onr)
>>> print(data.shape) # For the version 2 of PVAT dataset.
(7082, 22)
'''
(JD,ET,XLSE,YLSE,ZLSE,R,XME,YME,ZME,LAT,LON,ALT,XSCXL,XSCYL,XSCZL,YSCXL,YSCYL,YSCZL,SZA,ELV,NADIR,VELOC)=list(range(22))
#### In fact, (JD,ET,DATE,XLSE,YLSE.... ) is correct. 23 elements. However, numpy.loadtxt will not read DATE part.
#### So final product would be 22 elements as above.
def __init__(self, baseurl=None, rcfile=None):
if baseurl is None:
rc=Rc()
if rcfile is not None:
rc.append(rcfile)
baseurl=rc.get('cy1orb', 'pvatbase')
if baseurl is None:
baseurl='file:///Volumes/scidata/data/moon/saraopvat/'
self._baseurl=baseurl
logging.debug('PVAT instanced: base=%s'%baseurl)
self.clearCache()
[docs] def clearCache(self):
self._cache={}
[docs] def getPvatDataFrom(self, fname):
cols=list(range(23))
cols.pop(2)
# print cols
dat=numpy.loadtxt(fname, usecols=cols) # This reads yyyy-mm-ddThh:mm:ss incorrectly, so skipped.
logging.debug('File retrieved : %s'%str(dat.shape))
return dat
[docs] def getPvatOfOrbit(self, orbit):
'''Get a numpy array of the position for the specified orbit.
'''
if orbit in self._cache:
logging.debug('Cache data for orbit %d found.'%orbit)
return self._cache[orbit].copy()
uri=self._baseurl+'/pvat_%04d.gz'%orbit
logging.debug('Cache not found. Retrieve data from %s.'%uri)
f_ret=urllib.request.urlretrieve(uri)
logging.debug('File retrieved : %s'%f_ret[0])
self._cache[orbit]=self.getPvatDataFrom(f_ret[0])
return self._cache[orbit].copy()
[docs]class Cy1Pvat(Cy1PvatRaw):
''' Chandrayaan-1 PVAT file.
The PVAT file is in two format. Human readable version (master) and
pickled version (cache). This class handles both the file format.
>>> pvat = Cy1Pvat()
>>> onr = 1280
>>> data = pvat.getPvatOfOrbit(onr)
>>> print(data.shape)
(7082, 22)
>>> data = pvat.getPvatOfOrbit(onr)
If you do not want to use cache file, you may use :class:`Cy1PvatRaw` class.
The location of the PVAT master file should be specified by ``.irfpyrc`` file
with entry of ``cy1orb/pvatbase``.
The PVAT cache file should be specified with an entry of ``cy1orb/pvatcachedir``.
Note that if the ``pvatbase`` is specified as remote directory, no cache is used.
To remove the cache, you just use shell with
.. code-block:: sh
% rm /path/to/cache/pvat_????.pickle.gz
'''
def __init__(self, baseurl=None, rcfile=None):
Cy1PvatRaw.__init__(self, baseurl=baseurl, rcfile=rcfile)
rc = Rc()
if rcfile is not None:
rc.append(rcfile)
self.cachepath = rc.get('cy1orb', 'pvatcachedir')
def _searchcache(self, orbit):
''' Returns (cachefilename, cache_used)
'''
# logging.getLogger().setLevel(logging.DEBUG)
if self.cachepath is None:
logging.warn('Cache directory not found. Set "cy1orb/pvatcachedir" in .irfpyrc')
return None, False
cachefile = os.path.join(self.cachepath, 'pvat_%04d.pickle.gz' % orbit)
logging.debug('Cachefile: %s' % cachefile)
if os.path.exists(cachefile):
### Compare generation tim.e
uri=self._baseurl+'/pvat_%04d.gz'%orbit
logging.debug('Cache not found. Retrieve data from %s.'%uri)
f_ret=urllib.request.urlretrieve(uri)[0]
ctime_master = os.stat(f_ret).st_ctime
ctime_cache = os.stat(cachefile).st_ctime
if ctime_master > ctime_cache:
return cachefile, False
else:
return cachefile, False
return cachefile, True
[docs] def getPvatOfOrbit(self, orbit):
# logging.getLogger().setLevel(logging.DEBUG)
### If the data is in the memory
if orbit in self._cache:
logging.debug('Cache data for orbit %d in memoery' % orbit)
return Cy1PvatRaw.getPvatOfOrbit(self, orbit)
### First, search cache file.
cachefn, cacheused = self._searchcache(orbit)
### Open the cache
try:
fp = gzip.open(cachefn)
datarr = pickle.load(fp)
logging.debug('Data retrieved from pickle. %s' % str(datarr.shape))
except:
datarr = Cy1PvatRaw.getPvatOfOrbit(self, orbit)
logging.info(datarr)
try:
fp2 = gzip.open(cachefn, 'w')
pickle.dump(datarr, fp2)
fp2.close()
except Exception as e:
logging.info('Failed to write to a new cache file.')
logging.debug(e)
return datarr
[docs]class Cy1Orbit:
''' Class of CY-1 orbit information.
USAGE:
>>> cy1=Cy1Orbit()
>>> pos=cy1.getPosOfOrb(1234) # retrieve 4xN data (time,x,y,z) of the orbit 1234.
It will produce (4, N) data array representing (time, x, y, z) as the first argument.
'''
def __init__(self, baseurl=None, rcfile=None):
''' Create an instance of Cy1Orbit. Inside Cy1Pvat instance is stored.
'''
if baseurl is None:
rc=Rc()
if rcfile is not None:
rc.append(rcfile)
baseurl=rc.get('cy1orb', 'pvatbase')
if baseurl is None:
baseurl='file:///Volumes/scidata/data/moon/saraopvat/'
self._pvat = Cy1Pvat(baseurl=baseurl)
[docs] def loadOrbit(self, orbnr):
''' Load the orbit. '''
self._pvat.getPvatOfOrbit(orbnr)
[docs] def getPosOfOrb(self, orbit):
'''Get the spacecraft position list in LSE coordinates for the specified orbit.
:returns: 4xN array of the data. Time information is a Julday in double-precision.
'''
pvat=self._pvat.getPvatOfOrbit(orbit)
jd=pvat[:,Cy1Pvat.JD].copy()
x=pvat[:,Cy1Pvat.XLSE].copy()
y=pvat[:,Cy1Pvat.YLSE].copy()
z=pvat[:,Cy1Pvat.ZLSE].copy()
logging.debug('JD shape =%s / First data=%f'%(str(jd.shape), jd[0]))
logging.debug('X shape =%s / First data=%f'%(str(x.shape), x[0]))
logging.debug('Y shape =%s / First data=%f'%(str(y.shape), y[0]))
logging.debug('Z shape =%s / First data=%f'%(str(z.shape), z[0]))
return numpy.array([jd,x,y,z])
[docs] def getDayEqOfOrb2(self, orbit):
''' Get the dayside equator crossing time.
:param orbit: Orbit number.
:returns: A tuple that have (crossing_time, x, y, z).
``crossing_time`` is ``datetime.datetime`` instance.
'''
pvat = self._pvat.getPvatOfOrbit(orbit)
t = pvat[:, Cy1Pvat.JD].copy()
x = pvat[:, Cy1Pvat.XLSE].copy()
y = pvat[:, Cy1Pvat.YLSE].copy()
z = pvat[:, Cy1Pvat.ZLSE].copy()
td = numpy.ma.masked_where(x<0, t)
zd = abs(numpy.ma.masked_where(x<0, z))
idx = zd.argmin()
teq, xeq, yeq, zeq = t[idx], x[idx], y[idx], z[idx]
if abs(zeq) > 50:
raise RuntimeError('No equator crossing for the orbit %d??'%orbit)
jdeq = Julday(teq).getDatetime()
return (jdeq, xeq, yeq, zeq)
[docs] def getDayEqOfOrb(self, orbit):
''' Get the [jd_float, x, y, z] of the equator crossing in the dayside. z should be zero (not exactly due to implementation) and x should be plus.
'''
pvat =self._pvat.getPvatOfOrbit(orbit)
t=pvat[:,Cy1Pvat.JD].copy()
x=pvat[:,Cy1Pvat.XLSE].copy()
y=pvat[:,Cy1Pvat.YLSE].copy()
z=pvat[:,Cy1Pvat.ZLSE].copy()
ndat=len(t)
t0=t[0]
x0=x[0]
y0=y[0]
z0=z[0]
for i in range(1, ndat):
t1,x1,y1,z1=t[i],x[i],y[i],z[i]
if x1>0:
if z0*z1<=0:
if abs(z0)<abs(z1):
return (t0,x0,y0,z0)
else:
return (t1,x1,y1,z1)
t0,x0,y0,z0=t1,x1,y1,z1
raise RuntimeError('No equator crossing for the orbit %d??'%orbit)
[docs] def getPosOfOrb2(self, orbit):
'''Get the spacecraft position list in LSE coordinates for the specified orbit.
The method returns the orbit data as JdSeries.
The functionality is the same as getPosOfOrb() method, but the
difference is the convenience of the data processing with compensation of
execution time.
Compared to the getPosOfOrb() method, it takes 20-30 times more time.
:returns: JdSeries of the data. Data is a instance of irfpy.util.vector3d.Vector3d class.
'''
pvat=self._pvat.getPvatOfOrbit(orbit)
data=JdSeries()
for i in range(len(pvat[:,Cy1Pvat.JD])):
data.add(Julday(pvat[i,Cy1Pvat.JD]), Vector3d(pvat[i,Cy1Pvat.XLSE], pvat[i, Cy1Pvat.YLSE], pvat[i,Cy1Pvat.ZLSE]))
return data
[docs] def getNadirAngleOfOrb(self, orbit):
''' Get the spacecraft nadir direction of the specified orbit.
The spacecraft nadir direction is the angle between X_sc and exact nadir.
+X_sc is yaw, so, in nominal cases, it should be zero.
'''
pvat=self._pvat.getPvatOfOrbit(orbit)
jd=pvat[:,Cy1Pvat.JD].copy()
nad=pvat[:,Cy1Pvat.NADIR].copy()
return numpy.array([jd,nad])
[docs] def getSzaOfOrb(self, orbit):
''' Get the spacecraft solar zenith angle.
'''
pvat=self._pvat.getPvatOfOrbit(orbit)
jd=pvat[:,Cy1Pvat.JD].copy()
sza=pvat[:,Cy1Pvat.SZA].copy()
return numpy.array([jd,sza])
[docs] def getPosMeOfOrb(self, orbit):
''' Get position of spacecraft in the ME coordinate system.
:returns: 4xN array of the data. Time information is a Julday in double-precision.
'''
pvat=self._pvat.getPvatOfOrbit(orbit)
jd=pvat[:,Cy1Pvat.JD].copy()
x=pvat[:,Cy1Pvat.XME].copy()
y=pvat[:,Cy1Pvat.YME].copy()
z=pvat[:,Cy1Pvat.ZME].copy()
return numpy.array([jd,x,y,z])
[docs] def getSelenographicOfOrb(self, orbit):
''' Get position of spacecraft in the [lon, lat, height] list.
:returns: 4xN array of the data. Time information is a Julday in double-precision.
'''
pvat=self._pvat.getPvatOfOrbit(orbit)
jd=pvat[:,Cy1Pvat.JD].copy()
x=pvat[:,Cy1Pvat.LON].copy()
y=pvat[:,Cy1Pvat.LAT].copy()
z=pvat[:,Cy1Pvat.ALT].copy()
return numpy.array([jd,x,y,z])
[docs] def getSelenographicOfOrb2(self, orbit):
pvat=self._pvat.getPvatOfOrbit(orbit)
data=JdSeries()
for dat in pvat:
jd=Julday(dat[Cy1Pvat.JD])
data.add(jd, [dat[Cy1Pvat.LON], dat[Cy1Pvat.LAT], dat[Cy1Pvat.ALT]])
return data
[docs] def clearCache(self):
''' Clear cached data in Cy1Pvat.
'''
self._pvat.clearCache()
[docs]class Cy1Attitude:
''' Class of CY-1 attitude information.
USAGE::
> cy1att=Cy1Attitude()
> list = cy1att.getSc2LseMatrixOfOrb(942)
'''
def __init__(self, baseurl=None, rcfile=None):
''' Create an instance of Cy1Orbit. Inside Cy1Pvat instance is stored.
'''
if baseurl is None:
rc=Rc()
if rcfile is not None:
rc.append(rcfile)
baseurl=rc.get('cy1orb', 'pvatbase')
if baseurl is None:
baseurl='file:///Volumes/scidata/data/moon/saraopvat/'
self._pvat = Cy1Pvat(baseurl=baseurl)
[docs] def loadOrbit(self, orb):
self._pvat.getPvatOfOrbit(orb)
[docs] def getSc2LseMatrixOfOrb(self, orbit):
''' Get a list of the conversion matrix.
:returns: Matrix of SC to LSE conversion.
:rtype: ``JdSeries``
The conversion from SC to LSE is implemented.
Formulation is::
|Xlse| | XSCXL YSCXL ZSCXL | |Xsc|
|Ylse| = | XSCYL YSCYL ZSCYL | |Ysc|
|Zlse| | XSCZL YSCZL ZSCZL | |Zsc|
'''
pvat=self._pvat.getPvatOfOrbit(orbit)
list=JdSeries()
for i in range(len(pvat[:,0])):
jd=pvat[i,Cy1Pvat.JD]
xsxl=pvat[i,Cy1Pvat.XSCXL]
xsyl=pvat[i,Cy1Pvat.XSCYL]
xszl=pvat[i,Cy1Pvat.XSCZL]
ysxl=pvat[i,Cy1Pvat.YSCXL]
ysyl=pvat[i,Cy1Pvat.YSCYL]
yszl=pvat[i,Cy1Pvat.YSCZL]
xs=numpy.array([xsxl, xsyl, xszl])
ys=numpy.array([ysxl, ysyl, yszl])
zs=numpy.cross(xs,ys)
matx=numpy.array([xs,ys,zs]).transpose()
list.add(Julday(jd), matx)
return list
[docs] def clearCache(self):
''' Clear the cached data of Cy1Pvat
'''
self._pvat.clearCache()