First commit
Cees Bassa
parent
f854f7be9f
commit
b47057ad82
|
|
@ -0,0 +1,201 @@
|
|||
#!/usr/bin/env python
|
||||
import sys
|
||||
import numpy as np
|
||||
import cv2
|
||||
import time
|
||||
import ctypes
|
||||
import multiprocessing
|
||||
from astropy.time import Time
|
||||
from astropy.io import fits
|
||||
from sun import get_sunriseset
|
||||
|
||||
# Capture images
|
||||
def capture(buf,z1,t1,z2,t2,device,nx,ny,nz,tend):
|
||||
# Array flag
|
||||
first=True
|
||||
|
||||
# Loop until reaching end time
|
||||
while float(time.time())<tend:
|
||||
# Get frames
|
||||
for i in range(nz):
|
||||
# Get frame
|
||||
ret,frame=device.read()
|
||||
|
||||
# Skip lost frames
|
||||
if ret==True:
|
||||
# Store time
|
||||
t=float(time.time())
|
||||
|
||||
# Convert image to grayscale
|
||||
gray=np.asarray(cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)).astype(np.uint8)
|
||||
|
||||
# Store buffer
|
||||
z=gray
|
||||
|
||||
# Store results
|
||||
if first:
|
||||
z1[i]=z
|
||||
t1[i]=t
|
||||
else:
|
||||
z2[i]=z
|
||||
t2[i]=t
|
||||
|
||||
# Assign buffer ready
|
||||
if first:
|
||||
buf.value=1
|
||||
else:
|
||||
buf.value=2
|
||||
|
||||
# Swap flag
|
||||
first=not first
|
||||
|
||||
print("Exiting capture")
|
||||
|
||||
def compress(buf,z1,t1,z2,t2,nx,ny,nz,tend):
|
||||
# Flag to keep track of processed buffer
|
||||
process_buf=1
|
||||
|
||||
# Start processing
|
||||
while True:
|
||||
# Wait for buffer to become available
|
||||
while buf.value!=process_buf:
|
||||
time.sleep(1.0)
|
||||
|
||||
# Process first buffer
|
||||
if buf.value==1:
|
||||
t=t1
|
||||
z=z1.astype('float32')
|
||||
process_buf=2
|
||||
elif buf.value==2:
|
||||
t=t2
|
||||
z=z2.astype('float32')
|
||||
process_buf=1
|
||||
|
||||
# Compute statistics
|
||||
zmax=np.max(z,axis=0)
|
||||
znum=np.argmax(z,axis=0)
|
||||
zs1=np.sum(z,axis=0)-zmax
|
||||
zs2=np.sum(z*z,axis=0)-zmax*zmax
|
||||
zavg=zs1/float(nz-1)
|
||||
zstd=np.sqrt((zs2-zs1*zavg)/float(nz-2))
|
||||
|
||||
# Convert to float and flip
|
||||
zmax=np.flipud(zmax.astype('float32'))
|
||||
znum=np.flipud(znum.astype('float32'))
|
||||
zavg=np.flipud(zavg.astype('float32'))
|
||||
zstd=np.flipud(zstd.astype('float32'))
|
||||
|
||||
# Format time
|
||||
nfd="%s.%03d"%(time.strftime("%Y-%m-%dT%T", time.gmtime(t[0])),int((t[0]-np.floor(t[0]))*1000))
|
||||
t0=Time(nfd,format='isot')
|
||||
dt=t-t[0]
|
||||
|
||||
# Generate fits
|
||||
fname="%s.fits"%nfd
|
||||
|
||||
# Format header
|
||||
hdr=fits.Header()
|
||||
hdr['DATE-OBS']="%s"%nfd
|
||||
hdr['MJD-OBS']=t0.mjd
|
||||
hdr['EXPTIME']=dt[-1]-dt[0]
|
||||
hdr['NFRAMES']=nz
|
||||
hdr['CRPIX1']=float(nx)/2.0
|
||||
hdr['CRPIX2']=float(ny)/2.0
|
||||
hdr['CRVAL1']=0.0
|
||||
hdr['CRVAL2']=0.0
|
||||
hdr['CD1_1']=1.0/3600.0
|
||||
hdr['CD1_2']=0.0
|
||||
hdr['CD2_1']=0.0
|
||||
hdr['CD2_2']=1.0/3600.0
|
||||
hdr['CTYPE1']="RA---TAN"
|
||||
hdr['CTYPE2']="DEC--TAN"
|
||||
hdr['CUNIT1']="deg"
|
||||
hdr['CUNIT2']="deg"
|
||||
hdr['CRRES1']=0.0
|
||||
hdr['CRRES2']=0.0
|
||||
hdr['EQUINIX']=2000.0
|
||||
hdr['RADECSYS']="ICRS"
|
||||
hdr['COSPAR']=4171
|
||||
hdr['OBSERVER']="Cees Bassa"
|
||||
for i in range(nz):
|
||||
hdr['DT%04d'%i]=dt[i]
|
||||
for i in range(10):
|
||||
hdr['DUMY%03d'%i]=0.0
|
||||
|
||||
# Write fits file
|
||||
hdu=fits.PrimaryHDU(data=np.array([zavg,zstd,zmax,znum]),header=hdr)
|
||||
hdu.writeto(fname)
|
||||
print("Compressed",fname)
|
||||
|
||||
# Exit on end of capture
|
||||
if t[-1]>tend:
|
||||
break
|
||||
|
||||
# Exiting
|
||||
print("Exiting compress")
|
||||
|
||||
# Main function
|
||||
if __name__ == '__main__':
|
||||
# Current time
|
||||
tnow=float(time.time())
|
||||
tnow=1520573400.0
|
||||
tnow=1520510226.281413,
|
||||
|
||||
# Get sunrise and sunset times
|
||||
trise,tset=get_sunriseset(tnow,52.8344,6.3785,10,-6.0)
|
||||
|
||||
print(tnow,trise,tset)
|
||||
# Wait until sunset
|
||||
if (tset<trise) & (tnow<tset):
|
||||
print("Waiting for sunset at %s"%time.strftime("%FT%T",time.gmtime(tset)))
|
||||
|
||||
sys.exit()
|
||||
|
||||
|
||||
# Settings
|
||||
devid=int(sys.argv[1])
|
||||
nx=720
|
||||
ny=576
|
||||
nz=250
|
||||
|
||||
# Initialize device
|
||||
device=cv2.VideoCapture(devid)
|
||||
|
||||
# Set properties
|
||||
device.set(3,nx)
|
||||
device.set(4,ny)
|
||||
|
||||
# Initialize arrays
|
||||
z1base=multiprocessing.Array(ctypes.c_uint8,nx*ny*nz)
|
||||
z1=np.ctypeslib.as_array(z1base.get_obj()).reshape(nz,ny,nx)
|
||||
t1base=multiprocessing.Array(ctypes.c_double,nz)
|
||||
t1=np.ctypeslib.as_array(t1base.get_obj())
|
||||
z2base=multiprocessing.Array(ctypes.c_uint8,nx*ny*nz)
|
||||
z2=np.ctypeslib.as_array(z2base.get_obj()).reshape(nz,ny,nx)
|
||||
t2base=multiprocessing.Array(ctypes.c_double,nz)
|
||||
t2=np.ctypeslib.as_array(t2base.get_obj())
|
||||
buf=multiprocessing.Value('i',0)
|
||||
|
||||
|
||||
tend=float(time.time())+31.0
|
||||
|
||||
# Set processes
|
||||
pcapture=multiprocessing.Process(target=capture,args=(buf,z1,t1,z2,t2,device,nx,ny,nz,tend))
|
||||
pcompress=multiprocessing.Process(target=compress,args=(buf,z1,t1,z2,t2,nx,ny,nz,tend))
|
||||
|
||||
# Start
|
||||
pcapture.start()
|
||||
pcompress.start()
|
||||
|
||||
# End
|
||||
try:
|
||||
pcapture.join()
|
||||
pcompress.join()
|
||||
except KeyboardInterrupt:
|
||||
pcapture.terminate()
|
||||
pcompress.terminate()
|
||||
|
||||
# Release device
|
||||
device.release()
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,142 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import os
|
||||
from astropy.io import fits
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
from astropy import wcs
|
||||
from scipy import optimize
|
||||
|
||||
def residual(a,x,y,z):
|
||||
return z-(a[0]+a[1]*x+a[2]*y)
|
||||
|
||||
def offsets(ra0,de0,ra,de):
|
||||
w=wcs.WCS(naxis=2)
|
||||
w.wcs.crpix=np.array([0.0,0.0])
|
||||
w.wcs.crval=np.array([ra0,de0])
|
||||
w.wcs.cd=np.array([[1.0/3600.0,0.0],[0.0,1.0/3600.0]])
|
||||
w.wcs.ctype=["RA---TAN","DEC--TAN"]
|
||||
w.wcs.set_pv([(2,1,45.0)])
|
||||
|
||||
world=np.stack((ra,de),axis=-1)
|
||||
pix=w.wcs_world2pix(world,1)
|
||||
return pix[:,0],pix[:,1]
|
||||
|
||||
|
||||
def read_pixel_catalog(fname):
|
||||
d=np.loadtxt(fname)
|
||||
x,y,mag=d[:,0],d[:,1],d[:,2]
|
||||
|
||||
return x,y,mag
|
||||
|
||||
def read_tycho2_catalog(maxmag=9.0):
|
||||
hdu=fits.open("/home/bassa/code/python/satellite/astrometry-tycho2/build/tyc2.fits")
|
||||
|
||||
ra=hdu[1].data.field('RA')
|
||||
de=hdu[1].data.field('DEC')
|
||||
mag=hdu[1].data.field('MAG_VT')
|
||||
|
||||
c=mag<maxmag
|
||||
print("Selected %d stars from Tycho-2 catalog with m<%.1f"%(np.sum(c),maxmag))
|
||||
return ra[c],de[c],mag[c]
|
||||
|
||||
def match_catalogs(ra,de,x,y,xs,ys,rmin):
|
||||
res=[]
|
||||
for i in range(len(x)):
|
||||
dx,dy=xs-x[i],ys-y[i]
|
||||
r=np.sqrt(dx*dx+dy*dy)
|
||||
if np.min(r)<rmin:
|
||||
j=np.argmin(r)
|
||||
res.append([ra[i],de[i],xs[j],ys[j]])
|
||||
|
||||
return np.array(res)
|
||||
|
||||
def fit_wcs(res,x0,y0):
|
||||
ra,de,x,y=res[:,0],res[:,1],res[:,2],res[:,3]
|
||||
ra0,de0=np.mean(ra),np.mean(de)
|
||||
dx,dy=x-x0,y-y0
|
||||
|
||||
for i in range(5):
|
||||
w=wcs.WCS(naxis=2)
|
||||
w.wcs.crpix=np.array([0.0,0.0])
|
||||
w.wcs.crval=np.array([ra0,de0])
|
||||
w.wcs.cd=np.array([[1.0,0.0],[0.0,1.0]])
|
||||
w.wcs.ctype=["RA---TAN","DEC--TAN"]
|
||||
w.wcs.set_pv([(2,1,45.0)])
|
||||
|
||||
world=np.stack((ra,de),axis=-1)
|
||||
pix=w.wcs_world2pix(world,1)
|
||||
rx,ry=pix[:,0],pix[:,1]
|
||||
|
||||
ax,cov_q,infodict,mesg,ierr=optimize.leastsq(residual,[0.0,0.0,0.0],args=(dx,dy,rx),full_output=1)
|
||||
ay,cov_q,infodict,mesg,ierr=optimize.leastsq(residual,[0.0,0.0,0.0],args=(dx,dy,ry),full_output=1)
|
||||
|
||||
ra0,de0=w.wcs_pix2world([[ax[0],ay[0]]],1)[0]
|
||||
|
||||
w=wcs.WCS(naxis=2)
|
||||
w.wcs.crpix=np.array([x0,y0])
|
||||
w.wcs.crval=np.array([ra0,de0])
|
||||
w.wcs.cd=np.array([[ax[1],ax[2]],[ay[1],ay[2]]])
|
||||
w.wcs.ctype=["RA---TAN","DEC--TAN"]
|
||||
w.wcs.set_pv([(2,1,45.0)])
|
||||
|
||||
whdr={"CRPIX1":x0,"CRPIX2":y0,"CRVAL1":ra0,"CRVAL2":de0,"CD1_1":ax[1],"CD1_2":ax[2],"CD2_1":ay[1],"CD2_2":ay[2],"CTYPE1":"RA---TAN","CTYPE2":"DEC--TAN","CUNIT1":"DEG","CUNIT2":"DEG"}
|
||||
|
||||
return whdr
|
||||
|
||||
# FITS file
|
||||
fname="2018-02-12T17:39:55.270.fits"
|
||||
|
||||
# Read fits
|
||||
hdu=fits.open(fname)
|
||||
|
||||
# Fix to force WCS with two axes only
|
||||
hdu[0].header["NAXIS"]=2
|
||||
w=wcs.WCS(hdu[0].header)
|
||||
|
||||
# Read data
|
||||
data=hdu[0].data
|
||||
zavg,zstd,zmax,znum=data
|
||||
ny,nx=zavg.shape
|
||||
|
||||
# Get extrema
|
||||
vmin=np.mean(zavg)-2.0*np.std(zavg)
|
||||
vmax=np.mean(zavg)+3.0*np.std(zavg)
|
||||
|
||||
# Read pixel catalog
|
||||
xs,ys,ms=read_pixel_catalog(fname+".cat")
|
||||
|
||||
# Read Tycho2 catalog
|
||||
ra,de,mag=read_tycho2_catalog(9.0)
|
||||
|
||||
# Convert to pixel coordinates
|
||||
xc,yc=w.all_world2pix(ra,de,1)
|
||||
c=(xc>0) & (xc<nx) & (yc>0) & (yc<ny)
|
||||
xt,yt=xc[c],yc[c]
|
||||
rat,det=ra[c],de[c]
|
||||
|
||||
# Match catalogs
|
||||
res=match_catalogs(rat,det,xt,yt,xs,ys,10.0)
|
||||
|
||||
whdr=fit_wcs(res,nx//2,ny//2)
|
||||
|
||||
|
||||
|
||||
|
||||
#hdu=fits.PrimaryHDU(header=w.to_header(),data=data)
|
||||
#hdr=fits.Header()
|
||||
hdr=hdu[0].header
|
||||
for k,v in whdr.items():
|
||||
hdr[k]=v
|
||||
hdu=fits.PrimaryHDU(header=hdr,data=data)
|
||||
hdu.writeto("test.fits",overwrite=True)
|
||||
|
||||
# Convert catalog stars
|
||||
plt.figure(figsize=(10,8))
|
||||
plt.imshow(zavg,origin="lower",aspect=1.0,interpolation="None",vmin=vmin,vmax=vmax)
|
||||
plt.scatter(xs,ys,s=80,edgecolors="r",facecolors="none")
|
||||
plt.scatter(xt,yt,s=150,edgecolors="y",facecolors="none")
|
||||
#plt.scatter(xt[c],yt[c],s=50,edgecolors="k",facecolors="none")
|
||||
#plt.xlim(0,720)
|
||||
#plt.ylim(0,576)
|
||||
plt.show()
|
||||
|
|
@ -0,0 +1,366 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import sys,os,glob
|
||||
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+1e-9)
|
||||
|
||||
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()
|
||||
|
||||
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:])
|
||||
|
||||
# IOD position format 2: RA/DEC = HHMMmmm+DDMMmm MX (MX in minutes of arc)
|
||||
def format_position(ra,de):
|
||||
ram=60.0*ra/15.0
|
||||
rah=int(np.floor(ram/60.0))
|
||||
ram-=60.0*rah
|
||||
|
||||
des=np.sign(de)
|
||||
dem=60.0*np.abs(de)
|
||||
ded=int(np.floor(dem/60.0))
|
||||
dem-=60.0*ded
|
||||
|
||||
if des==-1:
|
||||
sign="-"
|
||||
else:
|
||||
sign="+"
|
||||
|
||||
return ("%02d%06.3f%c%02d%05.2f"%(rah,ram,sign,ded,dem)).replace(".","")
|
||||
|
||||
# Format IOD line
|
||||
def format_iod_line(norad,cospar,site_id,t,ra,de):
|
||||
pstr=format_position(ra,de)
|
||||
tstr=t.replace("-","").replace("T","").replace(":","").replace(".","")
|
||||
|
||||
return "%05d %-9s %04d G %s 17 25 %s 37 S"%(norad,cospar,site_id,tstr,pstr)
|
||||
|
||||
# Extract tracks
|
||||
def extract_tracks(fname,trkrmin,drdtmin,trksig,ntrkmin):
|
||||
# Read four frame
|
||||
ff=fourframe(fname)
|
||||
|
||||
# Skip saturated frames
|
||||
if np.sum(ff.zavg>240.0)/float(ff.nx*ff.ny)>0.95:
|
||||
return
|
||||
|
||||
# Read satelite IDs
|
||||
try:
|
||||
f=open(fname+".id")
|
||||
except OSError:
|
||||
print("Cannot open",fname+".id")
|
||||
else:
|
||||
lines=f.readlines()
|
||||
f.close()
|
||||
|
||||
# ppgplot arrays
|
||||
tr=np.array([-0.5,1.0,0.0,-0.5,0.0,1.0])
|
||||
heat_l=np.array([0.0,0.2,0.4,0.6,1.0])
|
||||
heat_r=np.array([0.0,0.5,1.0,1.0,1.0])
|
||||
heat_g=np.array([0.0,0.0,0.5,1.0,1.0])
|
||||
heat_b=np.array([0.0,0.0,0.0,0.3,1.0])
|
||||
|
||||
# Loop over identifications
|
||||
for line in lines:
|
||||
# Decode
|
||||
id=satid(line)
|
||||
|
||||
# Skip slow moving objects
|
||||
drdt=np.sqrt(id.dxdt**2+id.dydt**2)
|
||||
if drdt<drdtmin:
|
||||
continue
|
||||
|
||||
# Extract significant pixels
|
||||
x,y,t,sig=ff.significant(trksig,id.x0,id.y0,id.dxdt,id.dydt,trkrmin)
|
||||
|
||||
# Fit tracks
|
||||
if len(t)>ntrkmin:
|
||||
# Get times
|
||||
tmin=np.min(t)
|
||||
tmax=np.max(t)
|
||||
tmid=0.5*(tmax+tmin)
|
||||
mjd=ff.mjd+tmid/86400.0
|
||||
|
||||
# Skip if no variance in time
|
||||
if np.std(t-tmid)==0.0:
|
||||
continue
|
||||
|
||||
# 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)
|
||||
|
||||
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.pgopen("/xs")
|
||||
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))
|
||||
if (3600.0*ff.crres[0]<1e-3) | (3600.0*ff.crres[1]<1e-3) | (ff.crres[0]/ff.sx>2.0) | (ff.crres[1]/ff.sy>2.0):
|
||||
ppgplot.pgsci(2)
|
||||
else:
|
||||
ppgplot.pgsci(1)
|
||||
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.pgsci(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(ff.zmax,ff.nx,ff.ny,0,ff.nx-1,0,ff.ny-1,ff.vmax,ff.vmin,tr)
|
||||
ppgplot.pgbox("BCTSNI",0.,0,"BCTSNI",0.,0)
|
||||
ppgplot.pgstbg(1)
|
||||
|
||||
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([x0]),np.array([y0]),4)
|
||||
ppgplot.pgmove(xmin,ymin)
|
||||
ppgplot.pgdraw(xmax,ymax)
|
||||
ppgplot.pgsch(0.65)
|
||||
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
|
||||
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;
|
||||
|
||||
# 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__':
|
||||
# Minimum predicted velocity (pixels/s)
|
||||
drdtmin=10.0
|
||||
|
||||
# Track selection region around prediction (pixels)
|
||||
trkrmin=10.0
|
||||
|
||||
# Track selection sigma
|
||||
trksig=5.0
|
||||
|
||||
# Minimum track points
|
||||
ntrkmin=10
|
||||
|
||||
# extract_tracks("2018-02-26T05:26:15.801.fits",trkrmin,drdtmin,trksig,ntrkmin)
|
||||
|
||||
files=sorted(glob.glob("2*.fits"))
|
||||
|
||||
for file in files:
|
||||
extract_tracks(file,trkrmin,drdtmin,trksig,ntrkmin)
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,49 @@
|
|||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from astropy.io import fits
|
||||
import cv2
|
||||
import sys
|
||||
|
||||
hdu=fits.open(sys.argv[1])
|
||||
|
||||
data=hdu[0].data
|
||||
zavg,zstd,zmax,znum=data
|
||||
|
||||
zsig=(zmax-zavg)/(zstd+1e-9)
|
||||
|
||||
ny,nx=zsig.shape
|
||||
tmp=np.zeros(nx*ny*3,dtype='uint8').reshape(ny,nx,3)
|
||||
tmp[:,:,0]=zmax.astype('uint8')
|
||||
tmp[:,:,1]=0.1*zmax.astype('uint8')
|
||||
tmp[:,:,2]=0.1*zmax.astype('uint8')
|
||||
img=cv2.cvtColor(tmp,cv2.COLOR_BGR2RGB)
|
||||
|
||||
mask=np.array(zsig>5.0,dtype='uint8')
|
||||
tmp=np.zeros(nx*ny*3,dtype='uint8').reshape(ny,nx,3)
|
||||
tmp[:,:,0]=255*mask
|
||||
tmp[:,:,1]=255*mask
|
||||
tmp[:,:,2]=255*mask
|
||||
|
||||
mat=cv2.cvtColor(tmp,cv2.COLOR_BGR2GRAY)
|
||||
lines=cv2.HoughLines(mat,1,np.pi/180.0,200)
|
||||
if lines is not None:
|
||||
for rho,theta in lines[0]:
|
||||
a = np.cos(theta)
|
||||
b = np.sin(theta)
|
||||
x0 = a*rho
|
||||
y0 = b*rho
|
||||
x1 = int(x0 + 1000*(-b))
|
||||
y1 = int(y0 + 1000*(a))
|
||||
x2 = int(x0 - 1000*(-b))
|
||||
y2 = int(y0 - 1000*(a))
|
||||
|
||||
cv2.line(img,(x1,y1+10),(x2,y2+10),(0,0,255),2)
|
||||
|
||||
minLineLength = 100
|
||||
maxLineGap = 10
|
||||
lines = cv2.HoughLinesP(mat,1,np.pi/180,100,minLineLength,maxLineGap)
|
||||
if lines is not None:
|
||||
for x1,y1,x2,y2 in lines[0]:
|
||||
cv2.line(img,(x1,y1),(x2,y2),(0,255,0),2)
|
||||
|
||||
cv2.imwrite(sys.argv[1]+'.jpg',img)
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
class iod:
|
||||
"""IOD Observation"""
|
||||
|
||||
def __init__(self,line):
|
||||
s=line.split()
|
||||
self.iodline=line
|
||||
self.norad=int(s[0])
|
||||
self.cospar="%s %s"%(s[1],s[2])
|
||||
self.site=int(s[3])
|
||||
self.conditions=s[4]
|
||||
self.datestring=s[5]
|
||||
|
||||
|
||||
|
||||
def __repr__(self):
|
||||
return "%05d %s %04d %s %s"%(self.norad,self.cospar,self.site,self.conditions,self.datestring)
|
||||
|
||||
|
||||
iodlines=["10529 77 112D 4553 G 20120511214446298 17 25 2034286+472232 37 R",
|
||||
"10529 77 112D 4553 G 20120511214452938 17 25 2038653+454274 37 S",
|
||||
"23862 96 029D 4553 G 20120511221706217 17 25 2121225+480700 37 S",
|
||||
"23862 96 029D 4553 G 20120511221726241 17 25 2122590+453398 37 R"]
|
||||
|
||||
obslist=[iod(line) for line in iodlines]
|
||||
|
||||
print(type(obslist),type(obslist[0]))
|
||||
for obs in obslist:
|
||||
print(obs)
|
||||
|
|
@ -0,0 +1,69 @@
|
|||
#!/usr/bin/env python
|
||||
import numpy as np
|
||||
import sys
|
||||
from astropy.io import fits
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
|
||||
fname="2018-02-12T17:39:55.270.fits"
|
||||
|
||||
# Read fits
|
||||
hdu=fits.open(fname)
|
||||
|
||||
# Read data
|
||||
data=hdu[0].data
|
||||
zavg,zstd,zmax,znum=data
|
||||
|
||||
# Read header
|
||||
ny,nx=zavg.shape
|
||||
nz=hdu[0].header['NFRAMES']
|
||||
dt=[hdu[0].header['DT%04d'%i] for i in range(nz)]
|
||||
|
||||
# Sigma frame
|
||||
zsig=(zmax-zavg)/(zstd+1e-9)
|
||||
|
||||
# Position
|
||||
xm,ym=np.meshgrid(np.arange(nx),np.arange(ny))
|
||||
|
||||
# Selection width
|
||||
w=10
|
||||
sigma=5.0
|
||||
|
||||
# Selection
|
||||
c=(zsig>sigma) & (xm>w) & (xm<nx-w) & (ym>w) & (ym<ny-w)
|
||||
|
||||
x=np.ravel(xm[c])
|
||||
y=np.ravel(ym[c])
|
||||
inum=np.ravel(znum[c]).astype('int')
|
||||
sig=np.ravel(zsig[c])
|
||||
t=np.array([dt[i] for i in inum])
|
||||
|
||||
# Load predictions
|
||||
d=np.loadtxt(fname+".id",usecols=(1,2,3,4,5,6))
|
||||
|
||||
x0=d[:,0]
|
||||
y0=d[:,1]
|
||||
t0=np.zeros_like(x0)
|
||||
x1=d[:,2]
|
||||
y1=d[:,3]
|
||||
t1=d[:,4]
|
||||
|
||||
for i in range(len(x0)):
|
||||
xr=x0[i]+(x1[i]-x0[i])*(t-t0[i])/(t1[i]-t0[i])
|
||||
yr=y0[i]+(y1[i]-y0[i])*(t-t0[i])/(t1[i]-t0[i])
|
||||
r=np.sqrt((x-xr)**2+(y-yr)**2)
|
||||
c=r<10.0
|
||||
print(i,np.sum(c))
|
||||
plt.figure()
|
||||
plt.plot([x0[i],x1[i]],[y0[i],y1[i]],'r')
|
||||
plt.scatter(x[c],y[c],c=t[c],s=10.0*(sig[c]-sigma))
|
||||
plt.scatter(x,y,c=t,s=1)
|
||||
plt.show()
|
||||
|
||||
#fig=plt.figure()
|
||||
#ax=fig.add_subplot(111,projection='3d')
|
||||
#ax.scatter(x,y,t,c=t,s=10*(sig-5))
|
||||
|
||||
#for i in range(len(x0)):
|
||||
# ax.plot([x0[i],x1[i]],[y0[i],y1[i]],[t0[i],t1[i]],'k')
|
||||
#plt.show()
|
||||
|
|
@ -0,0 +1,191 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import numpy as np
|
||||
from astropy.io import fits
|
||||
from astropy.time import Time
|
||||
from astropy import wcs
|
||||
|
||||
class observation:
|
||||
"""Satellite observation"""
|
||||
|
||||
def __init__(self,ff,mjd,x0,y0):
|
||||
"""Define an observation"""
|
||||
|
||||
# Store
|
||||
self.mjd=mjd
|
||||
self.x0=x0
|
||||
self.y0=y0
|
||||
|
||||
# Get times
|
||||
self.nfd=Time(self.mjd,format='mjd',scale='utc').isot
|
||||
|
||||
# Correct for rotation
|
||||
tobs=Time(ff.mjd+0.5*ff.texp/86400.0,format='mjd',scale='utc')
|
||||
tobs.delta_ut1_utc=0
|
||||
hobs=tobs.sidereal_time("mean",longitude=0.0).degree
|
||||
tmid=Time(self.mjd,format='mjd',scale='utc')
|
||||
tmid.delta_ut1_utc=0
|
||||
hmid=tmid.sidereal_time("mean",longitude=0.0).degree
|
||||
|
||||
# Compute ra/dec
|
||||
world=ff.w.wcs_pix2world(np.array([[self.x0,self.y0]]),1)
|
||||
self.ra=world[0,0]+hobs-hmid
|
||||
self.de=world[0,1]
|
||||
|
||||
class satid:
|
||||
"""Satellite identifications"""
|
||||
|
||||
def __init__(self,line):
|
||||
s=line.split()
|
||||
self.nfd=s[0]
|
||||
self.x0=float(s[1])
|
||||
self.y0=float(s[2])
|
||||
self.t0=0.0
|
||||
self.x1=float(s[3])
|
||||
self.y1=float(s[4])
|
||||
self.t1=float(s[5])
|
||||
self.norad=int(s[6])
|
||||
self.catalog=s[7]
|
||||
self.state=s[8]
|
||||
self.dxdt=(self.x1-self.x0)/(self.t1-self.t0)
|
||||
self.dydt=(self.y1-self.y0)/(self.t1-self.t0)
|
||||
|
||||
def __repr__(self):
|
||||
return "%s %f %f %f -> %f %f %f %d %s %s"%(self.nfd,self.x0,self.y0,self.t0,self.x1,self.y1,self.t1,self.norad,self.catalog,self.state)
|
||||
|
||||
class fourframe:
|
||||
"""Four frame class"""
|
||||
|
||||
def __init__(self,fname=None):
|
||||
if fname==None:
|
||||
# Initialize empty fourframe
|
||||
self.nx=0
|
||||
self.ny=0
|
||||
self.nz=0
|
||||
self.mjd=-1
|
||||
self.nfd=None
|
||||
self.zavg=None
|
||||
self.zstd=None
|
||||
self.zmax=None
|
||||
self.znum=None
|
||||
self.dt=None
|
||||
self.site_id=0
|
||||
self.observer=None
|
||||
self.texp=0.0
|
||||
self.fname=None
|
||||
self.crpix=np.array([0.0,0.0])
|
||||
self.crval=np.array([0.0,0.0])
|
||||
self.cd=np.array([[1.0,0.0],[0.0,1.0]])
|
||||
self.ctype=["RA---TAN","DEC--TAN"]
|
||||
self.cunit=np.array(["deg","deg"])
|
||||
self.crres=np.array([0.0,0.0])
|
||||
else:
|
||||
# Read FITS file
|
||||
hdu=fits.open(fname)
|
||||
|
||||
# Read image planes
|
||||
self.zavg,self.zstd,self.zmax,self.znum=hdu[0].data
|
||||
|
||||
# Frame properties
|
||||
self.ny,self.nx=self.zavg.shape
|
||||
self.nz=hdu[0].header['NFRAMES']
|
||||
|
||||
# Read frame time oselfsets
|
||||
self.dt=np.array([hdu[0].header['DT%04d'%i] for i in range(self.nz)])
|
||||
|
||||
# Read header
|
||||
self.mjd=hdu[0].header['MJD-OBS']
|
||||
self.nfd=hdu[0].header['DATE-OBS']
|
||||
self.site_id=hdu[0].header['COSPAR']
|
||||
self.observer=hdu[0].header['OBSERVER']
|
||||
self.texp=hdu[0].header['EXPTIME']
|
||||
self.fname=fname
|
||||
|
||||
# Astrometry keywords
|
||||
self.crpix=np.array([hdu[0].header['CRPIX1'],hdu[0].header['CRPIX2']])
|
||||
self.crval=np.array([hdu[0].header['CRVAL1'],hdu[0].header['CRVAL2']])
|
||||
self.cd=np.array([[hdu[0].header['CD1_1'],hdu[0].header['CD1_2']],
|
||||
[hdu[0].header['CD2_1'],hdu[0].header['CD2_2']]])
|
||||
self.ctype=[hdu[0].header['CTYPE1'],hdu[0].header['CTYPE2']]
|
||||
self.cunit=[hdu[0].header['CUNIT1'],hdu[0].header['CUNIT2']]
|
||||
self.crres=np.array([hdu[0].header['CRRES1'],hdu[0].header['CRRES2']])
|
||||
|
||||
# Compute image properties
|
||||
self.sx=np.sqrt(self.cd[0,0]**2+self.cd[1,0]**2)
|
||||
self.sy=np.sqrt(self.cd[0,1]**2+self.cd[1,1]**2)
|
||||
self.wx=self.nx*self.sx
|
||||
self.wy=self.ny*self.sy
|
||||
self.vmin=np.mean(self.zmax)-2.0*np.std(self.zmax)
|
||||
self.vmax=np.mean(self.zmax)+6.0*np.std(self.zmax)
|
||||
|
||||
# Setup WCS
|
||||
self.w=wcs.WCS(naxis=2)
|
||||
self.w.wcs.crpix=self.crpix
|
||||
self.w.wcs.crval=self.crval
|
||||
self.w.wcs.cd=self.cd
|
||||
self.w.wcs.ctype=self.ctype
|
||||
self.w.wcs.set_pv([(2,1,45.0)])
|
||||
|
||||
|
||||
def significant(self,sigma,x0,y0,dxdt,dydt,rmin=10.0):
|
||||
"""Extract significant points"""
|
||||
|
||||
# Generate sigma frame
|
||||
zsig=(self.zmax-self.zavg)/(self.zstd+1e-9)
|
||||
|
||||
# Select
|
||||
c=(zsig>sigma)
|
||||
|
||||
# Positions
|
||||
xm,ym=np.meshgrid(np.arange(self.nx),np.arange(self.ny))
|
||||
x,y=np.ravel(xm[c]),np.ravel(ym[c])
|
||||
inum=np.ravel(self.znum[c]).astype('int')
|
||||
sig=np.ravel(zsig[c])
|
||||
t=np.array([self.dt[i] for i in inum])
|
||||
|
||||
# Predicted positions
|
||||
xr=x0+dxdt*t
|
||||
yr=y0+dydt*t
|
||||
r=np.sqrt((x-xr)**2+(y-yr)**2)
|
||||
c=r<rmin
|
||||
|
||||
return x[c],y[c],t[c],sig[c]
|
||||
|
||||
def track(self,dxdt,dydt,tref):
|
||||
"""Track and stack"""
|
||||
# Empty frame
|
||||
ztrk=np.zeros_like(self.zavg)
|
||||
|
||||
# Loop over frames
|
||||
for i in range(self.nz):
|
||||
dx=int(np.round(dxdt*(self.dt[i]-tref)))
|
||||
dy=int(np.round(dydt*(self.dt[i]-tref)))
|
||||
|
||||
# Skip if shift larger than image
|
||||
if np.abs(dx)>=self.nx:
|
||||
continue
|
||||
if np.abs(dy)>=self.ny:
|
||||
continue
|
||||
|
||||
# Extract range
|
||||
if dx>=0:
|
||||
i1min,i1max=dx,self.nx-1
|
||||
i2min,i2max=0,self.nx-dx-1
|
||||
else:
|
||||
i1min,i1max=0,self.nx+dx-1
|
||||
i2min,i2max=-dx,self.nx-1
|
||||
if dy>=0:
|
||||
j1min,j1max=dy,self.ny-1
|
||||
j2min,j2max=0,self.ny-dy-1
|
||||
else:
|
||||
j1min,j1max=0,self.ny+dy-1
|
||||
j2min,j2max=-dy,self.ny-1
|
||||
zsel=np.where(self.znum==i,self.zmax,0.0)
|
||||
ztrk[j2min:j2max,i2min:i2max]+=zsel[j1min:j1max,i1min:i1max]
|
||||
|
||||
return ztrk
|
||||
|
||||
def __repr__(self):
|
||||
return "%s %dx%dx%d %s %.3f %d %s"%(self.fname,self.nx,self.ny,self.nz,self.nfd,self.texp,self.site_id,self.observer)
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,45 @@
|
|||
#!/usr/bin/env python
|
||||
from astropy.coordinates import SkyCoord,EarthLocation,AltAz,get_sun
|
||||
from astropy.time import Time
|
||||
import astropy.units as u
|
||||
import numpy as np
|
||||
from scipy import optimize,interpolate
|
||||
import time
|
||||
|
||||
|
||||
def get_sunriseset(t,lat,lon,height,refalt):
|
||||
loc=EarthLocation(lat=lat*u.deg,lon=lon*u.deg,height=height*u.m)
|
||||
t0=Time(t,format='unix')
|
||||
print(t0.isot)
|
||||
# Compute solar position every 10 minutes
|
||||
t=Time(t0.mjd+np.arange(160)/160.0,format='mjd',scale='utc')
|
||||
psun=get_sun(t)
|
||||
hor=psun.transform_to(AltAz(obstime=t,location=loc))
|
||||
|
||||
# Interpolating function
|
||||
alt_diff=hor.alt.deg-refalt
|
||||
f=interpolate.interp1d(t.mjd,alt_diff)
|
||||
|
||||
# Find sunrise/sunset
|
||||
sign=alt_diff[1:]*alt_diff[:-1]
|
||||
idx=np.where(sign<0.0)[0]
|
||||
print(idx)
|
||||
for i in idx:
|
||||
# Set
|
||||
if alt_diff[i]>alt_diff[i+1]:
|
||||
tset=Time(optimize.bisect(f,t[i].mjd,t[i+1].mjd),format='mjd',scale='utc')
|
||||
# Rise
|
||||
else:
|
||||
trise=Time(optimize.bisect(f,t[i].mjd,t[i+1].mjd),format='mjd',scale='utc')
|
||||
|
||||
return trise.unix,tset.unix
|
||||
|
||||
if __name__ == '__main__':
|
||||
lat,lon,height=52.8344,6.3785,10
|
||||
|
||||
t=float(time.time())
|
||||
t=1520573400.0
|
||||
print(t,type(t))
|
||||
trise,tset=get_sunriseset(t,lat,lon,height,-6.0)
|
||||
|
||||
print(t,trise,tset)
|
||||
|
|
@ -0,0 +1,11 @@
|
|||
#!/usr/bin/env python
|
||||
from skyfield.positionlib import ICRF
|
||||
from skyfield.api import load,Topos
|
||||
|
||||
ts=load.timescale()
|
||||
t=ts.now()
|
||||
observer=Topos(latitude_degrees=52.8344,longitude_degrees=6.3785,elevation_m=10.0)
|
||||
p=ICRF([0.0,0.0,0.0],observer_data=observer,t=t)
|
||||
q=p.from_altaz(alt_degrees=30.0,az_degrees=30.0)
|
||||
|
||||
print(p,q)
|
||||
Loading…
Reference in New Issue