From 71940a9a4c00b677b7df25cb1b210132074231f8 Mon Sep 17 00:00:00 2001 From: Cees Bassa Date: Sat, 10 Mar 2018 14:09:16 +0100 Subject: [PATCH] Implemented acquisition on sunset to sunrise --- acquire.py | 49 +++++++++++++++++++---------- plot.py | 77 +++++++++++---------------------------------- utils.py | 92 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 144 insertions(+), 74 deletions(-) create mode 100644 utils.py diff --git a/acquire.py b/acquire.py index e86f136..d10b154 100644 --- a/acquire.py +++ b/acquire.py @@ -1,13 +1,16 @@ #!/usr/bin/env python +from __future__ import print_function import sys import numpy as np import cv2 import time import ctypes import multiprocessing +from astropy.coordinates import EarthLocation from astropy.time import Time from astropy.io import fits -from sun import get_sunriseset +import astropy.units as u +from utils import get_sunset_and_sunrise # Capture images def capture(buf,z1,t1,z2,t2,device,nx,ny,nz,tend): @@ -137,20 +140,37 @@ def compress(buf,z1,t1,z2,t2,nx,ny,nz,tend): # Main function if __name__ == '__main__': # Current time - tnow=float(time.time()) - tnow=1520573400.0 - tnow=1520510226.281413, + tnow=Time.now() + + # Set location + loc=EarthLocation(lat=52.8344*u.deg,lon=6.3785*u.deg,height=10*u.m) + + # Reference altitude + refalt=-6.0*u.deg # Get sunrise and sunset times - trise,tset=get_sunriseset(tnow,52.8344,6.3785,10,-6.0) + state,tset,trise=get_sunset_and_sunrise(tnow,loc,refalt) - print(tnow,trise,tset) - # Wait until sunset - if (tset=tset): + dt=np.floor((tset-tnow).to(u.s).value) + print("The sun is above the horizon. Waiting %.0f seconds."%dt) + tend=trise + try: + time.sleep(dt) + except KeyboardInterrupt: + sys.exit() + print("Starting data acquisition. Acquisition will end at "+tend.isot) # Settings devid=int(sys.argv[1]) @@ -176,12 +196,9 @@ if __name__ == '__main__': 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)) + pcapture=multiprocessing.Process(target=capture,args=(buf,z1,t1,z2,t2,device,nx,ny,nz,tend.unix)) + pcompress=multiprocessing.Process(target=compress,args=(buf,z1,t1,z2,t2,nx,ny,nz,tend.unix)) # Start pcapture.start() diff --git a/plot.py b/plot.py index d1118d7..bdc5a79 100644 --- a/plot.py +++ b/plot.py @@ -1,69 +1,30 @@ #!/usr/bin/env python + +import sys,os,glob +from stio import fourframe,satid,observation 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" +ff=fourframe("2018-02-26T03:09:05.866.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 +zsig=(ff.zmax-ff.zavg)/(ff.zstd+1e-9) sigma=5.0 -# Selection -c=(zsig>sigma) & (xm>w) & (xmw) & (ymsigma +print(np.sum(c)/float(ff.nx*ff.ny)) +xm,ym=np.meshgrid(np.arange(ff.nx),np.arange(ff.ny)) +x,y=np.ravel(xm[c]),np.ravel(ym[c]) +inum=np.ravel(ff.znum[c]).astype('int') sig=np.ravel(zsig[c]) -t=np.array([dt[i] for i in inum]) +t=np.array([ff.dt[i] for i in inum]) -# Load predictions -d=np.loadtxt(fname+".id",usecols=(1,2,3,4,5,6)) +plt.figure() +plt.scatter(x,y,c=t,s=1) +plt.colorbar() +plt.show() -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() +fig=plt.figure() +ax=fig.add_subplot(111,projection='3d') +ax.scatter(x,y,t,c=t,s=1) +plt.show() diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..d9fc679 --- /dev/null +++ b/utils.py @@ -0,0 +1,92 @@ +from astropy.coordinates import SkyCoord,EarthLocation,AltAz,FK5,get_sun +from astropy.time import Time +import astropy.units as u +import numpy as np +from scipy import interpolate + +# Sunrise/sunset algorithm from Astronomical Algorithms by Jean Meeus +def get_sunset_and_sunrise(tnow,loc,refalt): + # Get time + nmjd=64 + mjd0=np.floor(tnow.mjd) + mnow=tnow.mjd-mjd0 + mjd=np.linspace(mjd0-1.0,mjd0+2.0,nmjd) + t=Time(mjd,format='mjd',scale='utc') + + # Get sun position + psun=get_sun(t) + + # Correct for precession by converting to FK5 + pos=SkyCoord(ra=psun.ra.degree,dec=psun.dec.degree,frame='icrs',unit='deg').transform_to(FK5(equinox=t)) + + # Compute altitude extrema + de=np.mean(pos.dec) + minalt=np.arcsin(np.sin(loc.latitude)*np.sin(de)-np.cos(loc.latitude)*np.cos(de)) + maxalt=np.arcsin(np.sin(loc.latitude)*np.sin(de)+np.cos(loc.latitude)*np.cos(de)) + + # Never sets, never rises? + if minalt>refalt: + return "sun never sets",t[0],t[0] + elif maxalt