srg/srgcat/management/commands/00_init_skymaps.py
2024-04-26 12:43:00 +03:00

184 lines
6.0 KiB
Python

from django.core.management.base import BaseCommand, CommandError
from datetime import date
import datetime
from django.utils import timezone
import astropy
from astropy.io import ascii
import pandas as pd
import pymysql
from sqlalchemy import create_engine
from heasarc.tdat import tDat
from heasarc.models import HeasarcTable, TableColumn, HeasarcObjectClass, NSIDE_SOURCES, ORDER
from srgcat.models import ArtCat, ArtSource
from astropy.table import Table
from astropy_healpix import HEALPix, neighbours
from astropy.coordinates import SkyCoord # High-level coordinates
from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames
from astropy.coordinates import Angle, Latitude, Longitude # Angles
import astropy.units as u
from astropy.io import fits
from django.db.models import Q
from srgcat.models import SkyMaps
from monthplan.models import SurveyHealpixPlate, NSIDE_PLATES, ORDER
# XMMU JHHMMSS.s+/-DDMMSS.
# www.cosmos.esa.int/web/xmm-newton/source-naming-convention
def make_source_name(key, ra, dec):
c = SkyCoord(ra, dec, frame=FK5(), unit="deg")
str1 = c.to_string('hmsdms',alwayssign=False,pad=False,precision=1).split()
str2 = c.to_string('hmsdms',alwayssign=False,pad=False,precision=0).split()
name = key+" J%s%s" % (str1[0].replace('h','').replace('m','').replace('s',''),
str2[1].replace('d','').replace('m','').replace('s',''))
return name
def load_skymaps(filename):
#hdul = fits.open(filename)
#hdul.info()
table = Table.read(filename)
df = table.to_pandas()
for i in df.index:
skymap = SkyMaps(
SMAPNR = int(df['SMAPNR'][i]),
RA_MIN = float(df['RA_MIN'][i]),
RA_MAX = float(df['RA_MAX'][i]),
DE_MIN = float(df['DE_MIN'][i]),
DE_MAX = float(df['DE_MAX'][i]),
RA_CEN = float(df['RA_CEN'][i]),
DE_CEN = float(df['DE_CEN'][i]),
ELON_CEN = float(df['ELON_CEN'][i]),
ELAT_CEN = float(df['ELAT_CEN'][i]),
GLON_CEN = float(df['GLON_CEN'][i]),
GLAT_CEN = float(df['GLAT_CEN'][i]),
X_MIN = float(df['X_MIN'][i]),
Y_MIN = float(df['Y_MIN'][i]),
OWNER = int(df['OWNER'][i]),
N_NBRS = int(df['N_NBRS'][i]),
FIELD1 = int(df['FIELD1'][i]),
FIELD2 = int(df['FIELD2'][i]),
FIELD3 = int(df['FIELD3'][i]),
FIELD4 = int(df['FIELD4'][i]),
FIELD5 = int(df['FIELD5'][i]),
FIELD6 = int(df['FIELD6'][i]),
FIELD7 = int(df['FIELD7'][i]),
FIELD8 = int(df['FIELD8'][i]),
FIELD9 = int(df['FIELD9'][i]))
skymap.save()
pass
def find_neighbours():
skymaps = SkyMaps.objects.all()
for skymap in skymaps:
skymap.neighbours.clear()
try:
if(skymap.FIELD1 > 0):
n1 = SkyMaps.objects.get(SMAPNR=skymap.FIELD1)
skymap.neighbours.add(n1)
if(skymap.FIELD2 > 0):
n2 = SkyMaps.objects.get(SMAPNR=skymap.FIELD2)
skymap.neighbours.add(n2)
if(skymap.FIELD3 > 0):
n3 = SkyMaps.objects.get(SMAPNR=skymap.FIELD3)
skymap.neighbours.add(n3)
if(skymap.FIELD4 > 0):
n4 = SkyMaps.objects.get(SMAPNR=skymap.FIELD4)
skymap.neighbours.add(n4)
if(skymap.FIELD5 > 0):
n5 = SkyMaps.objects.get(SMAPNR=skymap.FIELD5)
skymap.neighbours.add(n5)
if(skymap.FIELD6 > 0):
n6 = SkyMaps.objects.get(SMAPNR=skymap.FIELD6)
skymap.neighbours.add(n6)
if(skymap.FIELD7 > 0):
n7 = SkyMaps.objects.get(SMAPNR=skymap.FIELD7)
skymap.neighbours.add(n7)
if(skymap.FIELD8 > 0):
n8 = SkyMaps.objects.get(SMAPNR=skymap.FIELD8)
skymap.neighbours.add(n8)
if(skymap.FIELD9 > 0):
n9 = SkyMaps.objects.get(SMAPNR=skymap.FIELD9)
skymap.neighbours.add(n9)
except:
print(skymap," neighbour not found")
skymap.save()
print(skymap)
def find_plates():
skymaps = SkyMaps.objects.all()
for skymap in skymaps:
skymap.survey_healpix_plate.clear()
plates = SurveyHealpixPlate.objects.all()
for plate in plates:
ra = plate.ra
dec = plate.dec
try:
skymap = SkyMaps.objects.get(Q(RA_MIN__lte=ra) & Q(RA_MAX__gt=ra) & Q(DE_MIN__lte=dec) & Q(DE_MAX__gt=dec))
except Exception as e:
print("%s (%s)" % (e.message, type(e)))
return
skymap.survey_healpix_plate.add(plate)
skymap.save()
around = neighbours(plate.healpix, NSIDE_PLATES, order=ORDER)
for index in around:
if(index<0):
continue
try:
pl = SurveyHealpixPlate.objects.get(healpix=index)
except:
print("Healpix index {} not found, exit".format(index))
return
skymap.survey_healpix_plate.add(pl)
skymap.save()
class Command(BaseCommand):
help = 'Initiates data dase'
def handle(self, *args, **options):
#skymaps = SkyMaps.objects.all()
#skymaps.delete()
#load_skymaps('/data/erosita/SKYMAPS.fits')
#find_neighbours()
find_plates()
return
"""
some checks:
"""
skymaps = SkyMaps.objects.all()[:10]
for skymap in skymaps:
print(skymap, skymap.FIELD1,
skymap.FIELD2,
skymap.FIELD3,
skymap.FIELD4,
skymap.FIELD5,
skymap.FIELD6,
skymap.FIELD7,
skymap.FIELD8,
skymap.FIELD9)
n = skymap.survey_healpix_plate.all()
for nn in n:
print(nn)
self.stdout.write(self.style.SUCCESS('Done'))