srg/monthplan/utils.py
2024-04-26 12:43:00 +03:00

276 lines
10 KiB
Python

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 <krivonos@cosmos.ru>'
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 <krivonos@cosmos.ru>
"""
# /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)