In [53]:
import os
import datetime
import numpy
import scipy.signal
from astropy.io import fits
import matplotlib.pyplot as plt
import matplotlib.dates as md
%matplotlib inline
In [48]:
paths = ['/home/roman/mnt/server-space/storage/bolidozor/ZVPP/ZVPP-R6/snapshots/2017/09/']
times = numpy.ndarray((0,2))
start_time = datetime.datetime.now()
fits_browsed = 0
for path in paths:
for root, dirs, files in os.walk(path):
print("")
print(root, " ")
for name in files:
if name.endswith(("snap.fits")):
hdulist = fits.open(os.path.join(root, name))
DATE_ts = datetime.datetime.strptime(hdulist[1].header['DATE'], '%Y-%m-%dT%H:%M:%S').timestamp()*1000+2*60*60*1000
crval = hdulist[1].header['CRVAL2']
#print(DATE_ts, crval)
time = [DATE_ts - hdulist[1].header['CDELT2']* hdulist[1].header['NAXIS2'], crval]
times = numpy.vstack( [times, time] )
hdulist.close()
print("+", end='')
fits_browsed += 1
times.sort(axis=0)
print("")
print("===================================")
print(fits_browsed, "was successfully processed")
print("It takes", datetime.datetime.now()-start_time)
In [68]:
paths = ['/home/roman/mnt/server-space/storage/bolidozor/ZVPP/ZVPP-R6/snapshots/2017/09/03/',
'/home/roman/mnt/server-space/storage/bolidozor/ZVPP/ZVPP-R6/snapshots/2017/09/04/',
'/home/roman/mnt/server-space/storage/bolidozor/ZVPP/ZVPP-R6/snapshots/2017/09/05/']
times_ts = numpy.ndarray((0,2))
start_time = datetime.datetime.now()
fits_browsed = 0
for path in paths:
for root, dirs, files in os.walk(path):
print("")
print(root, " ")
for name in files:
if name.endswith(("snap.fits")):
try:
hdulist = fits.open(os.path.join(root, name))
sysdate = hdulist[1].header['SYSDATE1']
sysdate_beg = sysdate - hdulist[1].header['CDELT2']* hdulist[1].header['NAXIS2']
crval = hdulist[1].header['CRVAL2']
time = [sysdate_beg, crval]
times_ts = numpy.vstack( [times_ts, time] )
hdulist.close()
print("+", end='')
fits_browsed += 1
except Exception:
print("-", end='')
times_ts.sort(axis=0)
print("")
print("===================================")
print(fits_browsed, "was successfully processed")
print("It takes", datetime.datetime.now()-start_time)
In [80]:
plt.figure(figsize=(30, 20))
data=md.date2num([datetime.datetime.fromtimestamp(ts, datetime.timezone.utc) for ts in times[:,0]/1000])
data_ts=md.date2num([datetime.datetime.fromtimestamp(ts, datetime.timezone.utc) for ts in times_ts[:,0]/1000])
plt.xticks( rotation=25 )
ax=plt.gca()
xfmt = md.DateFormatter('%Y-%m-%d %H:%M:%S')
ax.xaxis.set_major_formatter(xfmt)
ax.set_title('Difference between DATE and CRVAL2 (radio-observer time of 1st .FITS row)')
ax.set_xlabel('datetime [UTC]')
ax.set_ylabel('time difference (DATE - CRVAL2) [s]')
plt.plot(data_ts, (times_ts[:,0]-times_ts[:,1])/1000.0, 'or')
plt.plot(data, (times[:,0]-times[:,1])/1000.0, 'xb')
plt.plot(data, scipy.signal.savgol_filter(times[:,0]-times[:,1],501, 3)/1000.0, 'w')
plt.show()
In [18]:
fits_path = '/home/roman/mnt/server-space/storage/bolidozor/ZVPP/ZVPP-R6/snapshots/2017/09/04/19/20170904192530311_ZVPP-R6_snap.fits'
print("")
hdulist = fits.open(fits_path)
sysdate = hdulist[1].header['SYSDATE1']
#sysdate_beg = sysdate - hdulist[1].header['CDELT2']* hdulist[1].header['NAXIS2']
DATE_ts = datetime.datetime.strptime(hdulist[1].header['DATE-OBS'], '%Y-%m-%dT%H:%M:%S').timestamp()*1000.0
crval = hdulist[1].header['CRVAL2']
hdulist.close()
time = (DATE_ts - crval)/1000.0
if time>0:
print("difference between times is", time, "s. (SYSDATE is ahead, radio-observer time is late)")
else:
print("difference between times is", time, "s. (CRVAL2 is ahead, radio-observer time is in the future :-) )")
In [ ]: