2018-03-09 13:08:47 -07:00
|
|
|
#!/usr/bin/env python
|
2018-03-10 06:09:16 -07:00
|
|
|
from __future__ import print_function
|
2018-04-16 06:43:04 -06:00
|
|
|
import sys
|
|
|
|
import os
|
2018-03-09 13:08:47 -07:00
|
|
|
import numpy as np
|
|
|
|
import cv2
|
|
|
|
import time
|
|
|
|
import ctypes
|
|
|
|
import multiprocessing
|
2018-03-10 06:09:16 -07:00
|
|
|
from astropy.coordinates import EarthLocation
|
2018-03-09 13:08:47 -07:00
|
|
|
from astropy.time import Time
|
|
|
|
from astropy.io import fits
|
2018-03-10 06:09:16 -07:00
|
|
|
import astropy.units as u
|
2018-04-21 04:14:15 -06:00
|
|
|
from utils import get_sunset_and_sunrise
|
2018-03-30 06:09:09 -06:00
|
|
|
import logging
|
2018-04-21 04:14:15 -06:00
|
|
|
import configparser
|
2018-04-16 06:43:04 -06:00
|
|
|
|
2018-03-09 13:08:47 -07:00
|
|
|
# Capture images
|
2018-04-16 06:43:04 -06:00
|
|
|
def capture(buf, z1, t1, z2, t2, device, nx, ny, nz, tend):
|
2018-03-09 13:08:47 -07:00
|
|
|
# Array flag
|
2018-04-16 06:43:04 -06:00
|
|
|
first = True
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Loop until reaching end time
|
2018-04-16 06:43:04 -06:00
|
|
|
while float(time.time()) < tend:
|
2018-03-09 13:08:47 -07:00
|
|
|
# Get frames
|
|
|
|
for i in range(nz):
|
|
|
|
# Get frame
|
2018-04-16 06:43:04 -06:00
|
|
|
ret, frame = device.read()
|
|
|
|
|
2018-03-09 13:08:47 -07:00
|
|
|
# Skip lost frames
|
2018-04-16 06:43:04 -06:00
|
|
|
if ret is True:
|
2018-03-09 13:08:47 -07:00
|
|
|
# Store time
|
2018-04-16 06:43:04 -06:00
|
|
|
t = float(time.time())
|
|
|
|
|
2018-03-09 13:08:47 -07:00
|
|
|
# Convert image to grayscale
|
2018-04-16 06:43:04 -06:00
|
|
|
gray = np.asarray(cv2.cvtColor(
|
|
|
|
frame, cv2.COLOR_BGR2GRAY)).astype(np.uint8)
|
|
|
|
|
|
|
|
# Store buffer
|
|
|
|
z = gray
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Store results
|
|
|
|
if first:
|
2018-04-16 06:43:04 -06:00
|
|
|
z1[i] = z
|
|
|
|
t1[i] = t
|
2018-03-09 13:08:47 -07:00
|
|
|
else:
|
2018-04-16 06:43:04 -06:00
|
|
|
z2[i] = z
|
|
|
|
t2[i] = t
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Assign buffer ready
|
|
|
|
if first:
|
2018-04-16 06:43:04 -06:00
|
|
|
buf.value = 1
|
2018-03-09 13:08:47 -07:00
|
|
|
else:
|
2018-04-16 06:43:04 -06:00
|
|
|
buf.value = 2
|
|
|
|
|
2018-03-09 13:08:47 -07:00
|
|
|
# Swap flag
|
2018-04-16 06:43:04 -06:00
|
|
|
first = not first
|
2018-03-09 13:08:47 -07:00
|
|
|
|
2018-04-16 06:43:04 -06:00
|
|
|
logging.info("Exiting capture")
|
|
|
|
|
|
|
|
|
|
|
|
def compress(buf, z1, t1, z2, t2, nx, ny, nz, tend, path):
|
2018-03-09 13:08:47 -07:00
|
|
|
# Flag to keep track of processed buffer
|
2018-04-16 06:43:04 -06:00
|
|
|
process_buf = 1
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Start processing
|
|
|
|
while True:
|
|
|
|
# Wait for buffer to become available
|
2018-04-16 06:43:04 -06:00
|
|
|
while buf.value != process_buf:
|
2018-03-09 13:08:47 -07:00
|
|
|
time.sleep(1.0)
|
|
|
|
|
|
|
|
# Process first buffer
|
2018-04-16 06:43:04 -06:00
|
|
|
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
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Compute statistics
|
2018-04-16 06:43:04 -06:00
|
|
|
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))
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Convert to float and flip
|
2018-04-16 06:43:04 -06:00
|
|
|
zmax = np.flipud(zmax.astype('float32'))
|
|
|
|
znum = np.flipud(znum.astype('float32'))
|
|
|
|
zavg = np.flipud(zavg.astype('float32'))
|
|
|
|
zstd = np.flipud(zstd.astype('float32'))
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Format time
|
2018-04-16 06:43:04 -06:00
|
|
|
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]
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Generate fits
|
2018-04-16 06:43:04 -06:00
|
|
|
fname = "%s.fits" % nfd
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Format header
|
2018-04-16 06:43:04 -06:00
|
|
|
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['EQUINOX'] = 2000.0
|
|
|
|
hdr['RADECSYS'] = "ICRS"
|
|
|
|
hdr['COSPAR'] = cfg.getint('Common', 'observer_cospar')
|
|
|
|
hdr['OBSERVER'] = cfg.get('Common', 'observer_name')
|
2018-03-09 13:08:47 -07:00
|
|
|
for i in range(nz):
|
2018-04-16 06:43:04 -06:00
|
|
|
hdr['DT%04d' % i] = dt[i]
|
2018-03-09 13:08:47 -07:00
|
|
|
for i in range(10):
|
2018-04-16 06:43:04 -06:00
|
|
|
hdr['DUMY%03d' % i] = 0.0
|
|
|
|
|
2018-03-09 13:08:47 -07:00
|
|
|
# Write fits file
|
2018-04-16 06:43:04 -06:00
|
|
|
hdu = fits.PrimaryHDU(data=np.array([zavg, zstd, zmax, znum]),
|
|
|
|
header=hdr)
|
2018-04-21 04:14:15 -06:00
|
|
|
hdu.writeto(os.path.join(path,fname))
|
2018-04-16 06:43:04 -06:00
|
|
|
logging.info("Compressed %s" % fname)
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Exit on end of capture
|
2018-04-16 06:43:04 -06:00
|
|
|
if t[-1] > tend:
|
2018-03-09 13:08:47 -07:00
|
|
|
break
|
|
|
|
|
|
|
|
# Exiting
|
2018-03-30 06:09:09 -06:00
|
|
|
logging.info("Exiting compress")
|
2018-04-16 06:43:04 -06:00
|
|
|
|
|
|
|
|
|
|
|
# Main function
|
2018-03-09 13:08:47 -07:00
|
|
|
if __name__ == '__main__':
|
2018-04-21 04:14:15 -06:00
|
|
|
# Read commandline options (TODO; move to argparse?)
|
|
|
|
cfgfile = sys.argv[1]
|
|
|
|
if sys.argv[2] == "test":
|
|
|
|
testing = True
|
|
|
|
else:
|
|
|
|
testing = False
|
|
|
|
|
|
|
|
# Confiration Parsing
|
|
|
|
cfg = configparser.ConfigParser(inline_comment_prefixes=('#', ';'))
|
|
|
|
cfg.read(cfgfile)
|
|
|
|
|
|
|
|
# Get device id
|
|
|
|
devid = cfg.getint('Camera', 'device_id')
|
|
|
|
|
2018-03-09 13:08:47 -07:00
|
|
|
# Current time
|
2018-04-16 06:43:04 -06:00
|
|
|
tnow = Time.now()
|
2018-03-09 13:08:47 -07:00
|
|
|
|
2018-03-30 06:09:09 -06:00
|
|
|
# Get obsid
|
2018-04-16 06:43:04 -06:00
|
|
|
obsid = time.strftime("%Y%m%d_%H%M%S", time.gmtime())+"_%d" % devid
|
2018-03-30 06:09:09 -06:00
|
|
|
|
|
|
|
# Generate directory
|
2018-04-21 04:14:15 -06:00
|
|
|
path = os.path.join(cfg.get('Common', 'observations_path'),obsid)
|
2018-03-30 06:09:09 -06:00
|
|
|
os.makedirs(path)
|
2018-04-16 06:43:04 -06:00
|
|
|
|
2018-03-30 06:09:09 -06:00
|
|
|
# Setup logging
|
2018-04-21 04:14:15 -06:00
|
|
|
logging.basicConfig(filename=os.path.join(path,"acquire.log"), level=logging.DEBUG)
|
2018-03-30 06:09:09 -06:00
|
|
|
|
2018-03-10 06:09:16 -07:00
|
|
|
# Set location
|
2018-04-16 06:43:04 -06:00
|
|
|
loc = EarthLocation(lat=cfg.getfloat('Common', 'observer_lat')*u.deg,
|
|
|
|
lon=cfg.getfloat('Common', 'observer_lon')*u.deg,
|
|
|
|
height=cfg.getfloat('Common', 'observer_el')*u.m)
|
2018-03-09 13:08:47 -07:00
|
|
|
|
2018-04-21 04:14:15 -06:00
|
|
|
if testing == False:
|
|
|
|
# Reference altitude
|
|
|
|
refalt = -6.0 * u.deg
|
|
|
|
|
|
|
|
# Get sunrise and sunset times
|
|
|
|
state, tset, trise = get_sunset_and_sunrise(tnow, loc, refalt)
|
|
|
|
|
|
|
|
# Start/end logic
|
|
|
|
if state == "sun never rises":
|
|
|
|
logging.info("The sun never rises. Exiting program.")
|
2018-03-30 07:59:22 -06:00
|
|
|
sys.exit()
|
2018-04-21 04:14:15 -06:00
|
|
|
elif state == "sun never sets":
|
|
|
|
logging.info("The sun never sets.")
|
|
|
|
tend = tnow+24*u.h
|
|
|
|
elif (trise < tset):
|
|
|
|
logging.info("The sun is below the horizon.")
|
|
|
|
tend = trise
|
|
|
|
elif (trise >= tset):
|
|
|
|
dt = np.floor((tset-tnow).to(u.s).value)
|
|
|
|
logging.info("The sun is above the horizon. Sunset at %s." % tset.isot)
|
|
|
|
logging.info("Waiting %.0f seconds." % dt)
|
|
|
|
tend = trise
|
|
|
|
try:
|
|
|
|
time.sleep(dt)
|
|
|
|
except KeyboardInterrupt:
|
|
|
|
sys.exit()
|
|
|
|
else:
|
|
|
|
tend=tnow+31.0*u.s
|
|
|
|
|
2018-03-30 07:59:22 -06:00
|
|
|
|
2018-03-30 07:52:47 -06:00
|
|
|
logging.info("Starting data acquisition.")
|
|
|
|
logging.info("Acquisition will end at "+tend.isot)
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Settings
|
2018-04-21 04:14:15 -06:00
|
|
|
nx = cfg.getint('Camera', 'camera_x')
|
|
|
|
ny = cfg.getint('Camera', 'camera_y')
|
|
|
|
nz = cfg.getint('Camera', 'camera_frames')
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Initialize device
|
2018-04-16 06:43:04 -06:00
|
|
|
device = cv2.VideoCapture(devid)
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Set properties
|
2018-04-16 06:43:04 -06:00
|
|
|
device.set(3, nx)
|
|
|
|
device.set(4, ny)
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Initialize arrays
|
2018-04-16 06:43:04 -06:00
|
|
|
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)
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Set processes
|
2018-04-16 06:43:04 -06:00
|
|
|
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, path))
|
2018-03-09 13:08:47 -07:00
|
|
|
|
|
|
|
# Start
|
|
|
|
pcapture.start()
|
|
|
|
pcompress.start()
|
|
|
|
|
|
|
|
# End
|
|
|
|
try:
|
|
|
|
pcapture.join()
|
|
|
|
pcompress.join()
|
|
|
|
except KeyboardInterrupt:
|
|
|
|
pcapture.terminate()
|
|
|
|
pcompress.terminate()
|
2018-04-16 06:43:04 -06:00
|
|
|
|
2018-03-09 13:08:47 -07:00
|
|
|
# Release device
|
|
|
|
device.release()
|