sara_oplotΒΆ
Orbit plotter for CY1
USAGE: sara_oplot.py [-dAgg] [-h] [-v] [-b baseuri] [-o output] orbit_nr
- OPTS: -dAgg No display of window. Should be placed as first argument
- -h
help
- -v
verbose
- -b basename
Specify baseuri of the orbit data.
- -o output
Specify the output file name. If not specified (default), X window output.
- -m
Moon plot enabled.
- -p basename
Specify basuuri for Moon orbit data.
#!/usr/bin/python
'''Orbit plotter for CY1
USAGE: sara_oplot.py [-dAgg] [-h] [-v] [-b baseuri] [-o output] orbit_nr
OPTS: -dAgg No display of window. Should be placed as first argument
-h help
-v verbose
-b basename Specify baseuri of the orbit data.
-o output Specify the output file name. If not specified (default), X window output.
-m Moon plot enabled.
-p basename Specify basuuri for Moon orbit data.
'''
import os
import sys
argv = sys.argv
if len(argv) > 1 and argv[1] == '-dAgg':
import matplotlib
print("Using Agg.", file=sys.stderr)
matplotlib.use('Agg')
argv.remove('-dAgg')
from irfpy.util.julday import Julday
from irfpy.util import utc
from pylab import *
from matplotlib.patches import Arc
from numpy.ma import masked_where
import getopt
import datetime
import dateutil.parser
import logging
import subprocess
import re
from irfpy.cy1orb.Cy1Orbit import Cy1Orbit
def orbit_plot(axis, timelist, xlist, ylist, symbol):
print(len(timelist), len(xlist), len(ylist))
if len(timelist) != len(xlist) or len(timelist) != len(ylist):
raise RuntimeError('Data provided has inconsistency in length.\ntimelist=%d, xlist=%d, ylist=%d: but all 3 must be the same'%(len(timelist), len(xlist), len(ylist)))
return axis.plot(xlist, ylist, symbol)
class DataAndPlot:
'''Container class for data and plot.
'''
_hol_ax="+0"
_ver_ax="+1"
_dep_ax="+2"
_axis=None
def __init__(self, horizon_ax, vertical_ax, depth_ax):
''' Constructor specifying
horizon_ax, vertical_ax and depth_ax must be a string such as +x or -z and so on.
horizon_ax specifies which axis (x,y, or z) should be taken as a horizontal axis in the plot.
Similar way as vertical_ax.
depth_ax is how to decide in front of the plot plane.
In theory, depth_ax can be calculated as it complets the right hand frame, but here one must
specify whole the three axis.
Sample: Usual plot, horizon=+x, vertical=+y, depth=+z: __init__('+x','+y','+z')
Seen from North: horizon=+y, vertical=-x, depth=+z: __init__('+y','-x','+z')
'''
usable_strings=['+x','-x','+y','-y','+z','-z']
if not horizon_ax in usable_strings:
raise RuntimeError('horizon_ax must be specified from: '+str(usable_strings))
self._hor_ax=horizon_ax.replace('x','0').replace('y','1').replace('z','2')
if not vertical_ax in usable_strings:
raise RuntimeError('vertical_ax must be specified from: '+str(usable_strings))
self._ver_ax=vertical_ax.replace('x','0').replace('y','1').replace('z','2')
if not depth_ax in usable_strings:
raise RuntimeError('depth_ax must be specified from: '+str(usable_strings))
self._dep_ax=depth_ax.replace('x','0').replace('y','1').replace('z','2')
def plot_to_fig(self, fig, pos, timelist, poslist, title='Orbit plot', xlabel=None, ylabel=None):
'''Plot to figure.
fig is the Figure isntance, pos is the argument to be passed to fig.add_subplot(pos),
and timelist and poslist should be ndarray.
'''
if len(timelist)!=len(poslist):
raise RuntimeError('Data provided has inconsistency in length\ntimelist=%d, poslist=%d: but both are the same and must have more than 2')
Rm=1738.64
hor_sign = 1
if self._hor_ax[0]=='-':
hor_sign=-1
hor_nr = int(self._hor_ax[1])
if xlabel == None:
xlabel=self._hor_ax[1].replace('0', 'x [km]').replace('1', 'y [km]').replace('2', 'z [km]')
ver_sign = 1
if self._ver_ax[0]=='-':
ver_sign=-1
ver_nr = int(self._ver_ax[1])
if ylabel == None:
ylabel=self._ver_ax[1].replace('0', 'x [km]').replace('1', 'y [km]').replace('2', 'z [km]')
logging.info('HOR %d %d'%(hor_sign, hor_nr))
logging.info('VER %d %d'%(ver_sign, ver_nr))
self._axis=fig.add_subplot(pos, aspect='equal', xlim=[-2500*hor_sign,2500*hor_sign], ylim=[-2500*ver_sign,2500*ver_sign])
sphere = Arc((0,0), 2*Rm, 2*Rm, fill=False)
self._axis.add_patch(sphere)
orbit_plot(self._axis, timelist, poslist[:,hor_nr], poslist[:,ver_nr], 'b--')
self._axis.set_title(title)
self._axis.set_xlabel(xlabel)
self._axis.set_ylabel(ylabel)
dep_sign = 1
if self._dep_ax[0]=='-':
dep_sign=-1
dep_nr=int(self._dep_ax[1])
logging.info('DEP %d %d'%(dep_sign, dep_nr))
invisible=[(dep_sign)*poslist[i,dep_nr]<0 and (poslist[i,hor_nr]**2+poslist[i,ver_nr]**2)<Rm**2 for i in range(len(poslist[:,dep_nr]))]
hor_data = masked_where(invisible , poslist[:,hor_nr])
ver_data = masked_where(invisible , poslist[:,ver_nr])
line=self._axis.plot(hor_data, ver_data, 'b-', linewidth=2)
self._axis.text(-2150*hor_sign, -2150*ver_sign, timelist[0].strftime('%FT%T')+"--"+timelist[-1].strftime('%T'))
return self._axis
def set_limit(self):
'''Set the x and y limit (range).
By unknown reason, (probably auto range is enabled somewhere), the limit is not as expected.
Invoke the method at the end to modify the limit.
'''
hor_sign = 1
if self._hor_ax[0]=='-':
hor_sign=-1
ver_sign = 1
if self._ver_ax[0]=='-':
ver_sign=-1
self._axis.set_xlim(-2200*hor_sign,2200*hor_sign)
self._axis.set_ylim(-2200*ver_sign,2200*ver_sign)
def get_axis(self):
return self._axis
def ticks_plot(self, timelist, poslist):
logging.info('---- TICKS PLOT ----')
hor_nr = int(self._hor_ax[1])
ver_nr = int(self._ver_ax[1])
xlist = poslist[:,hor_nr]
ylist = poslist[:,ver_nr]
ticks_plot(self._axis, timelist, xlist, ylist)
def ticks_plot(axis, timelist, xlist, ylist):
logging.info("%d %d %d"%( len(timelist), len(xlist), len(ylist)))
if len(timelist) != len(xlist) or len(timelist) != len(ylist):
raise RuntimeError('Data provided has inconsistency in length.\ntimelist=%d, xlist=%d, ylist=%d: but all 3 must be the same'%(len(timelist), len(xlist), len(ylist)))
axis.plot(xlist, ylist, 'ro')
for i in range(len(timelist)):
axis.text(xlist[i], ylist[i], timelist[i].time())
def plot_data(timelist, poslist, verbose=False, Rm=1738.64, Rmp=1735.66, ticks=None, output=None, posttitle=''):
if verbose:
logging.basicConfig(level=logging.DEBUG)
logging.debug('DEBUG mode entered')
else:
logging.basicConfig(level=logging.INFO)
logging.info('Non DEBUG mode entered')
# Orbit data
logging.debug('position list (poslist) is converted to the numpy array.')
plist=array(poslist)
logging.debug(plist.shape)
if len(timelist)!=len(poslist):
raise RuntimeError('Data provided has inconsistency in length\ntimelist=%d, poslist=%d: but both are the same and must have more than 2')
timelist_t,plist_t=None,None
# Ticks data
if ticks!=None:
if len(ticks)!=2:
logging.warning('Specified ticks should have 2 component list: time and position. No plot is made.')
timelist_t,plist_t=None,None
else:
timelist_t, poslist_t = ticks
if len(timelist_t)!=len(poslist_t):
logging.warning('Specified ticks data has inconsistency in length of each data.\nlen(time)=%d len(pos)=%d\nNo plot is made.'%(len(timelist_t), len(poslist_t)))
timelist_t,plist_t=None,None
else:
plist_t = array(poslist_t)
d=arange(0, 181)*math.pi/180. # angle in radian for half circle
x=cos(d)*Rm # X-coords for half circle
y=sin(d)*Rm # Y-coords for half circle
fig=figure(figsize=(10,8))
### North view
nview=DataAndPlot('+y','-x','+z')
plt2=nview.plot_to_fig(fig, 222, timelist, plist, title='North view '+posttitle, xlabel='y [km]', ylabel='x [km]')
# nightside
xy=array([x,-y]).transpose()
night = Polygon(xy, closed=True, fill=True, alpha=0.2, facecolor='k')
plt2.add_patch(night)
if timelist_t != None:
nview.ticks_plot(timelist_t, plist_t)
### Dawn view
dview=DataAndPlot('+x','+z','-y')
plt3=dview.plot_to_fig(fig, 223, timelist, plist, title='Dawn view '+posttitle)
#night side
xy=array([-y,x]).transpose()
night = Polygon(xy, closed=True, fill=True, alpha=0.2, facecolor='k')
plt3.add_patch(night)
if timelist_t != None:
dview.ticks_plot(timelist_t, plist_t)
### Sun view
sview=DataAndPlot('+y','+z','+x')
plt4=sview.plot_to_fig(fig, 224, timelist, plist, title='Sun view '+posttitle)
if timelist_t != None:
sview.ticks_plot(timelist_t, plist_t)
nview.set_limit()
dview.set_limit()
sview.set_limit()
return (fig, (nview,dview,sview))
def _reformat(arrdata, dt=None):
tlist=[]
plist=[]
for i in range(len(arrdata[0,:])):
jd=Julday(arrdata[0,i]+0.5*Julday.SECOND).getDatetime()
jd=datetime.datetime(jd.year, jd.month, jd.day, jd.hour, jd.minute, jd.second, 0)
tlist.append(jd)
plist.append([ arrdata[1,i], arrdata[2,i], arrdata[3,i] ])
dtlist=[]
dplist=[]
if dt!=None:
# Search for 10 minutes ticks
tt = tlist[0]
microsecond=0
second=0
minute=int(tt.minute/10)*10
logging.debug('tt=%s'%str(tt))
t0 = datetime.datetime(tt.year, tt.month, tt.day, tt.hour, minute, second, microsecond)
logging.debug('t0=%s'%str(t0))
# Last data in the orbit
te=tlist[-1]
logging.debug('te=%s'%str(te))
while t0 < te:
if t0 in tlist:
logging.debug('%s %s %s'%(t0, tlist.index(t0), plist[(tlist.index(t0))]))
dtlist.append(t0)
dplist.append(plist[(tlist.index(t0))])
else:
logging.debug('Not found...%s'%t0)
t0=t0+datetime.timedelta(0, dt)
tlist=dtlist
plist=dplist
return (tlist, plist)
def oplot(args):
try:
opts, arg = getopt.getopt(args, 'hvb:o:p:m')
except getopt.GetoptError:
usage()
raise
basedir=None
output=None
withMoon=False
moonBasedir=None
for o,oa in opts:
if o in ('-h'):
usage()
return
if o in ('-v'):
logging.basicConfig(level=logging.DEBUG)
if o in ('-b'):
logging.debug('Basedir specified.')
basedir=oa
if o in ('-o'):
logging.debug('Output specified.')
output=oa
if o in ('-p'):
logging.debug('Basedir for Moon orbit set to %s'%oa)
moonBasedir=oa
if o in ('-m'):
logging.debug('Moon plot in upper-left')
withMoon=True
logging.info('Basedir=%s'%basedir)
logging.info('Output=%s'%output)
if len(arg) != 1:
usage()
return
else:
orbnr=int(arg[0])
logging.info('OrbitNr=%d'%orbnr)
if basedir==None:
cy1orb=Cy1Orbit()
else:
cy1orb=Cy1Orbit(baseurl=basedir)
arrdata = cy1orb.getPosOfOrb(orbnr)
dt=1200
timelist,poslist = _reformat(arrdata)
timelist_t, poslist_t = _reformat(arrdata, dt=dt)
fig, axes = plot_data(timelist, poslist, ticks=[timelist_t, poslist_t], output=output, posttitle='@%04d'%orbnr)
# fig, axes = plot_data(timelist, poslist)
if withMoon:
options=[]
if moonBasedir!=None:
options=['-b', moonBasedir]
options.append(timelist[0].strftime('%FT%T'))
import moon_oplot
moon_oplot.oplot(options, p_figure=fig, p_subplot='221')
if output==None:
print('plt.show() here')
plt.show()
fig.canvas.draw()
else:
print('savefig here')
plt.savefig(output)
def usage(fd=sys.stderr, exit=False):
print(__doc__, file=fd)
if __name__ == '__main__':
if len(argv)==1:
usage()
sys.exit(-1)
oplot(argv[1:])