spacecruft
/
sattools
1
0
Fork 0
Code Releases Activity

Track and stack added

pull/13/merge
Cees Bassa 2018-03-03 17:44:09 +01:00
parent cbebb43231
commit ad18c625dd
2 changed files with 208 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

211
python/extract_tracks.py
View File

@ -5,16 +5,93 @@ from stio import fourframe,satid,observation
import numpy as np
import ppgplot
import matplotlib.pyplot as plt
from scipy import optimize,ndimage
from astropy import wcs
from astropy.coordinates import SkyCoord
# Gaussian model
def model(a,nx,ny):
x,y=np.meshgrid(np.arange(nx),np.arange(ny))
dx,dy=(x-a[0])/a[2],(y-a[1])/a[2]
arg=-0.5*(dx**2+dy**2)
return a[3]*np.exp(arg)+a[4]
# Residual function
def residual(a,img):
ny,nx=img.shape
mod=model(a,nx,ny)
return (img-mod).ravel()
# Find peak
def peakfind(img,w=1.0):
# Find approximate location
ny,nx=img.shape
i=np.argmax(img)
y0=int(i/nx)
x0=i-y0*nx
# Image properties
imgavg=np.mean(img)
imgstd=np.std(img)
# Estimate
a=np.array([x0,y0,w,img[y0,x0]-imgavg,imgavg])
q,cov_q,infodict,mesg,ier=optimize.leastsq(residual,a,args=(img),full_output=1)
# Extract
xc,yc,w=q[0],q[1],q[2]
# Significance
sigma=(a[3]-imgavg)/imgstd
return xc,yc,w,sigma
# Plot selection
def plot_selection(id,x0,y0,dt=2.0,w=10.0):
dx,dy=id.x1-id.x0,id.y1-id.y0
ang=np.arctan2(dy,dx)
r=np.sqrt(dx**2+dy**2)
drdt=r/(id.t1-id.t0)
sa,ca=np.sin(ang),np.cos(ang)
dx=np.array([-dt,-dt,dt,dt,-dt])*drdt
dy=np.array([w,-w,-w,w,w])
x=ca*dx-sa*dy+x0
y=sa*dx+ca*dy+y0
ppgplot.pgsci(7)
ppgplot.pgline(x,y)
return
# Check if point is inside selection
def inside_selection(id,x0,y0,xmid,ymid,dt=2.0,w=10.0):
dx,dy=id.x1-id.x0,id.y1-id.y0
ang=-np.arctan2(dy,dx)
r=np.sqrt(dx**2+dy**2)
drdt=r/(id.t1-id.t0)
sa,ca=np.sin(ang),np.cos(ang)
dx,dy=x0-xmid,y0-ymid
rm=ca*dx-sa*dy
wm=sa*dx+ca*dy
dtm=rm/drdt
if (abs(wm)<w) & (abs(dtm)<dt):
return True
else:
return False
# Get COSPAR ID
def get_cospar(norad):
f=open(os.getenv("ST_DATADIR")+"/data/desig.txt")
lines=f.readlines()
f.close()
cospar=([line for line in lines if str(norad) in line])[0].split()[1]
try:
cospar=([line for line in lines if str(norad) in line])[0].split()[1]
except IndexError:
cospar="18500A"
return "%2s %s"%(cospar[0:2],cospar[2:])
@ -79,8 +156,26 @@ def extract_tracks(fname,trkrmin,drdtmin,trksig,ntrkmin):
# Fit tracks
if len(t)>ntrkmin:
obs=observation(ff,x,y,t,sig)
# Get times
tmin=np.min(t)
tmax=np.max(t)
tmid=0.5*(tmax+tmin)
mjd=ff.mjd+tmid/86400.0
# Very simple polynomial fit; no weighting, no cleaning
px=np.polyfit(t-tmid,x,1)
py=np.polyfit(t-tmid,y,1)
# Extract results
x0,y0=px[1],py[1]
dxdt,dydt=px[0],py[0]
xmin=x0+dxdt*(tmin-tmid)
ymin=y0+dydt*(tmin-tmid)
xmax=x0+dxdt*(tmax-tmid)
ymax=y0+dydt*(tmax-tmid)
cospar=get_cospar(id.norad)
obs=observation(ff,mjd,x0,y0)
iod_line="%s"%format_iod_line(id.norad,cospar,ff.site_id,obs.nfd,obs.ra,obs.de)
print(iod_line)
@ -98,6 +193,7 @@ def extract_tracks(fname,trkrmin,drdtmin,trksig,ntrkmin):
# Plot
ppgplot.pgopen(fname.replace(".fits","")+"_%05d.png/png"%id.norad)
#ppgplot.pgopen("/xs")
ppgplot.pgpap(0.0,1.0)
ppgplot.pgsvp(0.1,0.95,0.1,0.8)
@ -122,25 +218,116 @@ def extract_tracks(fname,trkrmin,drdtmin,trksig,ntrkmin):
ppgplot.pgsci(4)
elif id.catalog.find("inttles.tle")>0:
ppgplot.pgsci(3)
ppgplot.pgpt(np.array([obs.x0]),np.array([obs.y0]),4)
ppgplot.pgmove(obs.xmin,obs.ymin)
ppgplot.pgdraw(obs.xmax,obs.ymax)
ppgplot.pgpt(np.array([x0]),np.array([y0]),4)
ppgplot.pgmove(xmin,ymin)
ppgplot.pgdraw(xmax,ymax)
ppgplot.pgsch(0.65)
ppgplot.pgtext(np.array([obs.x0]),np.array([obs.y0])," %05d"%id.norad)
ppgplot.pgtext(np.array([x0]),np.array([y0])," %05d"%id.norad)
ppgplot.pgsch(1.0)
ppgplot.pgsci(1)
ppgplot.pgend()
elif id.catalog.find("classfd.tle")>0:
# Track and stack
# else:
# ztrk=ff.track(id.dxdt,id.dydt,0.0)
t=np.linspace(0.0,ff.texp)
x,y=id.x0+id.dxdt*t,id.y0+id.dydt*t
c=(x>0) & (x<ff.nx) & (y>0) & (y<ff.ny)
# Skip if no points selected
if np.sum(c)==0:
continue
# Compute track
tmid=np.mean(t[c])
mjd=ff.mjd+tmid/86400.0
xmid=id.x0+id.dxdt*tmid
ymid=id.y0+id.dydt*tmid
ztrk=ndimage.gaussian_filter(ff.track(id.dxdt,id.dydt,tmid),1.0)
vmin=np.mean(ztrk)-2.0*np.std(ztrk)
vmax=np.mean(ztrk)+6.0*np.std(ztrk)
# Select region
xmin=int(xmid-100)
xmax=int(xmid+100)
ymin=int(ymid-100)
ymax=int(ymid+100)
if xmin<0: xmin=0
if ymin<0: ymin=0
if xmax>ff.nx: xmax=ff.nx-1
if ymax>ff.ny: ymax=ff.ny-1
# Find peak
x0,y0,w,sigma=peakfind(ztrk[ymin:ymax,xmin:xmax])
x0+=xmin
y0+=ymin
# Skip if peak is not significant
if sigma<trksig:
continue
# Skip if point is outside selection area
if inside_selection(id,xmid,ymid,x0,y0)==False:
continue;
# plt.figure(figsize=(15,10))
# plt.imshow(ztrk,origin='lower')
# plt.scatter(id.x0,id.y0,s=150,edgecolors="y",facecolors="none")
# plt.show()
# Format IOD line
cospar=get_cospar(id.norad)
obs=observation(ff,mjd,x0,y0)
iod_line="%s"%format_iod_line(id.norad,cospar,ff.site_id,obs.nfd,obs.ra,obs.de)
print(iod_line)
if id.catalog.find("classfd.tle")>0:
outfname="classfd.dat"
elif id.catalog.find("inttles.tle")>0:
outfname="inttles.dat"
else:
outfname="catalog.dat"
f=open(outfname,"a")
f.write("%s\n"%iod_line);
f.close()
# Plot
ppgplot.pgopen(fname.replace(".fits","")+"_%05d.png/png"%id.norad)
ppgplot.pgpap(0.0,1.0)
ppgplot.pgsvp(0.1,0.95,0.1,0.8)
ppgplot.pgsch(0.8)
ppgplot.pgmtxt("T",6.0,0.0,0.0,"UT Date: %.23s COSPAR ID: %04d"%(ff.nfd,ff.site_id))
ppgplot.pgmtxt("T",4.8,0.0,0.0,"R.A.: %10.5f (%4.1f'') Decl.: %10.5f (%4.1f'')"%(ff.crval[0],3600.0*ff.crres[0],ff.crval[1],3600.0*ff.crres[1]))
ppgplot.pgmtxt("T",3.6,0.0,0.0,"FoV: %.2f\\(2218)x%.2f\\(2218) Scale: %.2f''x%.2f'' pix\\u-1\\d"%(ff.wx,ff.wy,3600.0*ff.sx,3600.0*ff.sy))
ppgplot.pgmtxt("T",2.4,0.0,0.0,"Stat: %5.1f+-%.1f (%.1f-%.1f)"%(np.mean(ff.zmax),np.std(ff.zmax),ff.vmin,ff.vmax))
ppgplot.pgmtxt("T",0.3,0.0,0.0,iod_line)
ppgplot.pgsch(1.0)
ppgplot.pgwnad(0.0,ff.nx,0.0,ff.ny)
ppgplot.pglab("x (pix)","y (pix)"," ")
ppgplot.pgctab(heat_l,heat_r,heat_g,heat_b,5,1.0,0.5)
ppgplot.pgimag(ztrk,ff.nx,ff.ny,0,ff.nx-1,0,ff.ny-1,vmax,vmin,tr)
ppgplot.pgbox("BCTSNI",0.,0,"BCTSNI",0.,0)
ppgplot.pgstbg(1)
plot_selection(id,xmid,ymid)
ppgplot.pgsci(0)
if id.catalog.find("classfd.tle")>0:
ppgplot.pgsci(4)
elif id.catalog.find("inttles.tle")>0:
ppgplot.pgsci(3)
ppgplot.pgpt(np.array([id.x0]),np.array([id.y0]),17)
ppgplot.pgmove(id.x0,id.y0)
ppgplot.pgdraw(id.x1,id.y1)
ppgplot.pgpt(np.array([x0]),np.array([y0]),4)
ppgplot.pgsch(0.65)
ppgplot.pgtext(np.array([id.x0]),np.array([id.y0])," %05d"%id.norad)
ppgplot.pgsch(1.0)
ppgplot.pgsci(1)
ppgplot.pgend()
if __name__ == '__main__':

31
python/stio.py
View File

@ -8,30 +8,17 @@ from astropy import wcs
class observation:
"""Satellite observation"""
def __init__(self,ff,x,y,t,sig):
"""Fit a satellite track"""
# Get times
self.tmin=np.min(t)
self.tmax=np.max(t)
self.tmid=0.5*(self.tmin+self.tmax)
self.mjd=ff.mjd+self.tmid/86400.0
self.nfd=Time(self.mjd,format='mjd',scale='utc').isot
# Very simple polynomial fit; no weighting, no cleaning
px=np.polyfit(t-self.tmid,x,1)
py=np.polyfit(t-self.tmid,y,1)
def __init__(self,ff,mjd,x0,y0):
"""Define an observation"""
# Store
self.x0=px[1]
self.y0=py[1]
self.dxdt=px[0]
self.dydt=py[0]
self.xmin=self.x0+self.dxdt*(self.tmin-self.tmid)
self.ymin=self.y0+self.dydt*(self.tmin-self.tmid)
self.xmax=self.x0+self.dxdt*(self.tmax-self.tmid)
self.ymax=self.y0+self.dydt*(self.tmax-self.tmid)
self.mjd=mjd
self.x0=x0
self.y0=y0
# Get times
self.nfd=Time(self.mjd,format='mjd',scale='utc').isot
# Correct for rotation
hobs=Time(ff.mjd+0.5*ff.texp/86400.0,format='mjd',scale='utc').sidereal_time("mean",longitude=0.0).degree
hmid=Time(self.mjd,format='mjd',scale='utc').sidereal_time("mean",longitude=0.0).degree