import numpy as np from math import pi, cos, sin from Quaternion import Quat from dateutil.parser import parse from datetime import timedelta import glob import re import fnmatch import os import logging from scipy.spatial.transform import Rotation, Slerp import pandas as pd from pandas import ExcelWriter from pandas import ExcelFile from astropy_healpix import HEALPix from astropy import units as u from astropy.coordinates import SkyCoord, RangeError from astropy.coordinates import ICRS, Galactic, FK5 from astropy.coordinates import Angle import astropy.units as u from astropy.time.formats import erfa from astropy.time import Time, TimeDelta, TimezoneInfo, TimeFromEpoch from astropy.io import fits from astropy.table import Table from django.contrib.auth.models import User from django.core import mail from django.template.loader import render_to_string from django.utils.html import strip_tags from django.core.mail import send_mail from monthplan.models import Survey, SurveyPath, NSIDE_PLATES, ORDER, SurveyHealpixPlate from monthplan.models import Head, Observation, Seance, Correction, Scan from srglib.utils import date2mission, slerp, str2date, TZ_MSK #MJDREF = 51543.875 #TZ_UTC = TimezoneInfo(utc_offset=0*u.hour) #TZ_MSK = TimezoneInfo(utc_offset=3*u.hour) def logic(index,nskip): if index % nskip == 0: return False return True def load_surveypath_data(file, nskip=100, forced=False, notify=False): match = '20??????_??????_???????????.iki' filename=os.path.basename(file) if fnmatch.fnmatch(filename, match): logging.debug("Match filename: Read {}".format(file)) else: logging.debug("Match filename: Skip {}".format(file)) obsid = re.findall("\d+", file)[2] if(obsid == '11000400200'): """ Merge two files for 11000400200: cat /mnt/nfs/npol/orientation/DMV/ARJ-SCAN/20200218_220000_11000300500.iki /mnt/nfs/npol/orientation/DMV/ARJ-SCAN/20200219_180000_11000400200.iki > /srv/srg-plan/srg/data/ARJ-SCAN/20200219_180000_11000400200.actual.iki """ file='/export/django/srg/data/npol/ARJ-SCAN/20200219_180000_11000400200.actual.iki' logging.debug("Special case for 11000400200, load correct file {}".format(file)) try: survey=Survey.objects.get(experiment__exact=obsid) except: logging.debug("Survey {} is not found".format(obsid)) return """ Update eroday start/stop """ survey.eroday_start=date2mission(survey.start)/14400 survey.eroday_stop=date2mission(survey.stop)/14400 survey.mjd_start = Time(survey.start , format='datetime').mjd survey.mjd_stop = Time(survey.stop , format='datetime').mjd survey.save() if(forced == True): survey.loaded=False survey.save() if(survey.loaded == True): logging.debug("Survey {} is already loaded, skip.".format(obsid)) return logging.info("Load {}, for Survey {}".format(file,obsid)) spath=survey.surveypath_set.all() spath.delete() csvfile=file df = pd.read_csv(csvfile,names=['date', 'time', 'q1', 'q2','q3','q4','dummy'], header=None, delim_whitespace=True, skipfooter=1, engine='python', skiprows = lambda x: logic(x, nskip)) for i in df.index: dtstr="%s %s" % (df['date'][i], df['time'][i]) dt = str2date(dtstr, tzone=TZ_MSK) ts = date2mission(dt) eroday=ts/14400 dtime = Time(dt , format='datetime') mjd=dtime.mjd q1=df['q1'][i] q2=df['q2'][i] q3=df['q3'][i] q4=df['q4'][i] quat=Quat(attitude=[q2,q3,q4,q1]) sp=SurveyPath(survey=survey,dtime=dt,mjd=mjd,obt=ts,eroday=eroday,q1=q1,q2=q2,q3=q3,q4=q4,ra=quat.ra,dec=quat.dec,roll=quat.roll) sp.save() survey.loaded=True """ Mark survey as loaded """ survey.save() tolerance=10 hp = HEALPix(nside=NSIDE_PLATES, order=ORDER, frame=FK5()) surveypaths = survey.surveypath_set.all() for i in range(len(surveypaths)-1): delta=(surveypaths[i+1].dtime-surveypaths[i].dtime).total_seconds() if(abs(int(delta)) > tolerance): nstep = int(abs(delta) / tolerance) delta_arr = np.linspace(0,1,nstep) quat = slerp([surveypaths[i].q2, surveypaths[i].q3, surveypaths[i].q4, surveypaths[i].q1], [surveypaths[i+1].q2, surveypaths[i+1].q3, surveypaths[i+1].q4, surveypaths[i+1].q1], delta_arr) for q in quat: qfin=Quat(attitude=q) crd = SkyCoord(qfin.ra, qfin.dec, frame="fk5", unit="deg") heal = hp.skycoord_to_healpix(crd) try: plate = SurveyHealpixPlate.objects.get(healpix=heal) except: logging.error('Error: SurveyHealpixPlate not found, run ./manage.py init_survey_healpix_plates first') break queryset = survey.surveyhealpixplate_set.all() if not queryset.filter(pk=plate.pk).exists(): plate.survey.add(survey) plate.save() if(notify == True): users = User.objects.filter(groups__name='srg-iki-admin') to_emails=[] for user in users: to_emails.append(user.email) subject = "SRG-ArXiv: New SRG flight plan {}".format(obsid) html_message = render_to_string('monthplan/surveypath_loaded_email.html', {'survey':survey,}) plain_message = strip_tags(html_message) from_email = 'Roman Krivonos ' for email in to_emails: mail.send_mail(subject, plain_message, from_email, [email], html_message=html_message) def load_monthplan_fits(filename): table = Table.read(filename) df = table.to_pandas() #print("Column headings:") #print(df.columns.str.encode('utf-8')) hdul = fits.open(filename) #hdul.info() version=int(hdul[1].header['VERSION']) """ delete existing month plan at this month """ dt = parse(hdul[1].header['START']) middle=dt+timedelta(days=10) heads = Head.objects.all() for h in heads: dt = h.get_datetime() if(dt.year == middle.year and dt.month == middle.month): print("Found {} v{} (new v{})".format(h,h.version,version)) if(h.version < version): print("Delete {} v{} (new v{})".format(h,h.version,version)) h.delete() else: print("Skip loading, already exists: {} v{}".format(h,h.version)) return print("Loading ",filename) head=Head() head.start=hdul[1].header['START'] head.stop=hdul[1].header['STOP'] head.gentime=hdul[1].header['GENTIME'] head.version=hdul[1].header['VERSION'] head.author=hdul[1].header['AUTHOR'] head.file=os.path.basename(filename) head.nrows=df.shape[0] head.save() for i in df.index: row=(i+1) if(df['TYPE'][i].decode().strip() == 'SEANCE'): """ print(df['START'][i].decode().strip()) """ seance=Seance(head=head, row=row) seance.start=df['START'][i].decode().strip() seance.stop=df['STOP'][i].decode().strip() seance.stations=df['STATIONS'][i].decode().strip() seance.guid=df['GUID'][i].decode().strip() seance.save() if(df['TYPE'][i].decode().strip() == 'OBSERVATION'): obs=Observation(head=head, row=row) obs.start=df['START'][i].decode().strip() obs.stop=df['STOP'][i].decode().strip() obs.target=df['TARGET'][i].decode().strip() obs.experiment=df['EXPERIMENT'][i].decode().strip() obs.ra=df['RA'][i] obs.dec=df['DEC'][i] obs.ra=df['RA'][i] obs.roll_angle=df['ROLL_ANGLE'][i] obs.sun_xoz_angle=df['SUN_XOZ_ANGLE'][i] obs.guid=df['GUID'][i].decode().strip() obs.save() if(df['TYPE'][i].decode().strip() == 'CORRECTION'): cor=Correction(head=head, row=row) cor.start=df['START'][i].decode().strip() cor.stop=df['STOP'][i].decode().strip() cor.impstart=df['IMPSTART'][i].decode().strip() cor.guid=df['GUID'][i].decode().strip() cor.save() if(df['TYPE'][i].decode().strip() == 'SCAN'): scan=Scan(head=head, row=row) scan.start=df['START'][i].decode().strip() scan.stop=df['STOP'][i].decode().strip() scan.target=df['TARGET'][i].decode().strip() scan.experiment=df['EXPERIMENT'][i].decode().strip() scan.ra=df['RA'][i] scan.dec=df['DEC'][i] scan.roll_angle=df['ROLL_ANGLE'][i] scan.sun_xoz_angle=df['SUN_XOZ_ANGLE'][i] scan.template=int(df['TEMPLATE'][i].decode().strip()) scan.guid=df['GUID'][i].decode().strip() scan.save() if(df['TYPE'][i].decode().strip() == 'SURVEY'): sur=Survey(head=head, row=row) sur.start=df['START'][i].decode().strip() sur.stop=df['STOP'][i].decode().strip() sur.target=df['TARGET'][i].decode().strip() sur.experiment=df['EXPERIMENT'][i].decode().strip() sur.ra_p=df['RA_P'][i] sur.dec_p=df['DEC_P'][i] sur.ra_z0=df['RA_Z0'][i] sur.dec_z0=df['DEC_Z0'][i] sur.ra_zk=df['RA_ZK'][i] sur.dec_zk=df['DEC_ZK'][i] sur.z_speed=df['Z_SPEED'][i] sur.guid=df['GUID'][i].decode().strip() sur.save() def load_monthplan_dir(dir): """ Loads monthplan from dir Parameters ---------- dir : str Absolute path to the input directory :Author: Roman Krivonos """ # /export/django/srg/data/npol/fits/RG_MonthPlan_2020-06_v07.fits logging.info("Loading monthplan from {}".format(dir)) files = glob.glob("{}/RG_MonthPlan_202?-??_v??.fits".format(dir)) files.sort(key=os.path.getmtime) for filename in files: #print(filename) load_monthplan_fits(filename)