New python-based data acquisition script
Cees Bassa
parent
c02d9d26e2
commit
f33a2ec958
|
|
@ -0,0 +1,176 @@
|
|||
#!/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
|
||||
|
||||
# Capture images
|
||||
def capture(buf,z1,t1,z2,t2,device,nx,ny,nz):
|
||||
# Array flag
|
||||
first=True
|
||||
|
||||
# For ever loop
|
||||
while True:
|
||||
# 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
|
||||
|
||||
def compress(buf,z1,t1,z2,t2,nx,ny,nz):
|
||||
# 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)
|
||||
|
||||
# Main function
|
||||
if __name__ == '__main__':
|
||||
# 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)
|
||||
|
||||
pcapture=multiprocessing.Process(target=capture,args=(buf,z1,t1,z2,t2,device,nx,ny,nz))
|
||||
pcompress=multiprocessing.Process(target=compress,args=(buf,z1,t1,z2,t2,nx,ny,nz))
|
||||
|
||||
pcapture.start()
|
||||
pcompress.start()
|
||||
pcapture.join()
|
||||
pcompress.join()
|
||||
|
||||
# print("Starting acquisition")
|
||||
# for i in range(5):
|
||||
# Capture
|
||||
# tstart=time.time()
|
||||
# t,z=capture(device,nx,ny,nz)
|
||||
# print("Captured %d %dx%d frames in %.3f s"%(nz,nx,ny,time.time()-tstart))
|
||||
|
||||
# Compress
|
||||
# tstart=time.time()
|
||||
# compress(t,z)
|
||||
# print("Compressed %d %dx%d frames in %.3f s"%(nz,nx,ny,time.time()-tstart))
|
||||
|
||||
# Release device
|
||||
device.release()
|
||||
|
||||
|
||||
|
|
@ -1,75 +0,0 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import numpy as np
|
||||
import cv2
|
||||
import time
|
||||
import subprocess
|
||||
import os
|
||||
|
||||
# Make directory
|
||||
if not os.path.exists("/dev/shm/video0"):
|
||||
os.mkdir("/dev/shm/video0")
|
||||
|
||||
# Get sunset and sunrise times (sattools allnight program)
|
||||
tsolar=subprocess.check_output(['allnight']).replace("\n","").split(" ")
|
||||
|
||||
# Convert to time structs
|
||||
tset=time.strptime(tsolar[0]+" UTC", "%Y-%m-%dT%H:%M:%S %Z")
|
||||
trise=time.strptime(tsolar[1]+" UTC", "%Y-%m-%dT%H:%M:%S %Z")
|
||||
tnow=time.gmtime()
|
||||
dtset=time.mktime(tset)-time.mktime(tnow)
|
||||
dtrise=time.mktime(trise)-time.mktime(tnow)
|
||||
|
||||
# Wait for sunset
|
||||
#if dtset>0:
|
||||
# print time.strftime("%FT%T",time.gmtime())+" waiting for sunset (%ds)"%dtset
|
||||
# time.sleep(dtset)
|
||||
|
||||
# Settings
|
||||
device=cv2.VideoCapture(0)
|
||||
device.set(3,720)
|
||||
device.set(4,576)
|
||||
|
||||
# Set counter
|
||||
iframe=1
|
||||
|
||||
# Start capture
|
||||
print time.strftime("%FT%T",time.gmtime())+" start capture"
|
||||
while dtrise>0:
|
||||
# Get frame
|
||||
ret,frame=device.read()
|
||||
|
||||
# Skip lost frames
|
||||
if ret==True:
|
||||
# Get time
|
||||
t=float(time.time())
|
||||
|
||||
# Format time
|
||||
nfd="%s.%03d"%(time.strftime("%Y-%m-%dT%T", time.gmtime(t)),int((t-np.floor(t))*1000))
|
||||
|
||||
# Get Size
|
||||
ny,nx=frame.shape[0],frame.shape[1]
|
||||
|
||||
# Convert image to grayscale
|
||||
gray=np.asarray(cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)).astype(np.int8)
|
||||
|
||||
# Open output file
|
||||
f=open("/dev/shm/video0/img%06d.pgm"%iframe,"w")
|
||||
f.write("P5\n# %s\n%d %d\n255\n"%(nfd,nx,ny))
|
||||
f.write(gray)
|
||||
f.close()
|
||||
|
||||
# Log
|
||||
if iframe%250==1:
|
||||
print "%06d: %s %dx%d"%(iframe,nfd,nx,ny)
|
||||
|
||||
# Increment
|
||||
iframe+=1
|
||||
|
||||
# Refresh time
|
||||
tnow=time.gmtime()
|
||||
dtrise=time.mktime(trise)-time.mktime(tnow)
|
||||
|
||||
|
||||
# End capture
|
||||
print time.strftime("%FT%T",time.gmtime())+" end capture"
|
||||
|
|
@ -1,135 +0,0 @@
|
|||
#!/usr/bin/env python
|
||||
import numpy as np
|
||||
from astropy.time import Time
|
||||
from astropy.io import fits
|
||||
import time
|
||||
import glob
|
||||
import os
|
||||
|
||||
def read_pgm(fname):
|
||||
# Open file
|
||||
f=open(fname,"r")
|
||||
|
||||
# Read lines
|
||||
line1=f.readline()
|
||||
line2=f.readline()
|
||||
line3=f.readline()
|
||||
line4=f.readline()
|
||||
|
||||
# Read parameters
|
||||
nfd=line2.split(" ")[1]
|
||||
nx,ny=int(line3.split(" ")[0]),int(line3.split(" ")[1])
|
||||
|
||||
# Read image
|
||||
z=np.fromfile(f,dtype='uint8',count=nx*ny).astype('float32')
|
||||
|
||||
# Close file
|
||||
f.close()
|
||||
|
||||
return z,nfd
|
||||
|
||||
# Create fits file
|
||||
def create_fits_file(nx,ny,nz,path,iframe):
|
||||
# Allocate
|
||||
z=np.empty(nz*ny*nx,dtype='float32').reshape(nz,nx*ny)
|
||||
mjd=np.empty(nz,dtype='float64')
|
||||
dt=np.empty(nz,dtype='float32')
|
||||
|
||||
# Loop over frames
|
||||
tstart=time.time()
|
||||
for i in xrange(nz):
|
||||
z[i],nfd=read_pgm(path+"/img%06d.pgm"%(iframe+i))
|
||||
|
||||
# Format time
|
||||
t=Time(nfd,format='isot')
|
||||
if i==0:
|
||||
t0=t
|
||||
mjd[i]=t.mjd
|
||||
dt[i]=86400.0*(mjd[i]-mjd[0])
|
||||
print "Files read in %.3f s"%(time.time()-tstart)
|
||||
|
||||
# Compute statistics
|
||||
tstart=time.time()
|
||||
zmax=np.max(z,axis=0)
|
||||
znum=np.argmax(z,axis=0)
|
||||
z1=np.sum(z,axis=0)-zmax
|
||||
z2=np.sum(z*z,axis=0)-zmax*zmax
|
||||
zavg=z1/float(nz-1)
|
||||
zstd=np.sqrt((z2-z1*zavg)/float(nz-2))
|
||||
print "Statistics in %.3f s"%(time.time()-tstart)
|
||||
|
||||
# Reshape, reformat and flip
|
||||
tstart=time.time()
|
||||
zmax=np.flipud(zmax.astype('float32').reshape(ny,nx))
|
||||
znum=np.flipud(znum.astype('float32').reshape(ny,nx))
|
||||
zavg=np.flipud(zavg.astype('float32').reshape(ny,nx))
|
||||
zstd=np.flipud(zstd.astype('float32').reshape(ny,nx))
|
||||
z=np.array([zavg,zstd,zmax,znum])
|
||||
print "Reshape in %.3f s"%(time.time()-tstart)
|
||||
|
||||
# Filename
|
||||
fname="%s.fits"%t0
|
||||
|
||||
# Format header
|
||||
hdr=fits.Header()
|
||||
hdr['DATE-OBS']="%s"%t0
|
||||
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
|
||||
hdr['CD1_2']=0.0
|
||||
hdr['CD2_1']=0.0
|
||||
hdr['CD2_2']=1.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 xrange(nz):
|
||||
hdr['DT%04d'%i]=dt[i]
|
||||
for i in xrange(10):
|
||||
hdr['DUMY%03d'%i]=0.0
|
||||
|
||||
# Write fits file
|
||||
hdu=fits.PrimaryHDU(data=z,header=hdr)
|
||||
hdu.writeto(fname,clobber=True)
|
||||
|
||||
return fname
|
||||
|
||||
# Main function
|
||||
if __name__ == '__main__':
|
||||
# Settings
|
||||
nx=720
|
||||
ny=576
|
||||
nz=250
|
||||
path="/dev/shm/video0"
|
||||
|
||||
# Start forever loop
|
||||
while True:
|
||||
# Find files
|
||||
files=sorted(glob.glob(path+"/img??????.pgm"))
|
||||
|
||||
# Enough files
|
||||
if len(files)>nz:
|
||||
# Get first frame
|
||||
iframe=int(files[0].replace(path+"/img","").replace(".pgm",""))
|
||||
|
||||
tstart=time.time()
|
||||
fname=create_fits_file(nx,ny,nz,path,iframe)
|
||||
tend=time.time()
|
||||
print "Created %s in %.2f"%(fname,tend-tstart)
|
||||
|
||||
# Remove files
|
||||
for i in xrange(nz):
|
||||
os.remove(path+"/img%06d.pgm"%(iframe+i))
|
||||
else:
|
||||
time.sleep(1)
|
||||
Loading…
Reference in New Issue