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

684 lines
23 KiB
Python

import random
import math
from django.db.models import DecimalField, Value
from django.db.models import IntegerField, BooleanField
from django.db.models import Max, Subquery, OuterRef, Exists, Min, Case, When, Min, Func
from math import isfinite, sqrt
from PIL import Image
from io import BytesIO
import pylab
import numpy as np
from astroquery.skyview import SkyView
import requests
import logging
import astropy
from astropy.io import fits
from astropy.io import ascii
from astropy_healpix import HEALPix
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_healpix import HEALPix, neighbours
from astropy import units as au
from astropy.time.formats import erfa
from astropy.time import Time, TimeDelta, TimezoneInfo, TimeFromEpoch
from astropy.table import Table
from django.db.models import Q
from django.db.models import FloatField, Value
from artsurvey.models import MetaSource
from artsurvey.models import ArtSurvey, EnergyBand
from artsurvey.models import ArtSurveySource
from artsurvey.models import ArtSurveyParams
from erosurvey.models import NSIDE_PLATES, ORDER_PLATES
from heasarc.models import NSIDE_SOURCES, ORDER
from heasarc.models import HeasarcBase
from heasarc.models import HeasarcSwiftBAT105m
from heasarc.models import HeasarcIntRefCat
from heasarc.models import HeasarcXTEASSCAT
from heasarc.models import HeasarcXMMSL2
from heasarc.models import Heasarc4FGL
from heasarc.models import HeasarcIntegral2020
from heasarc.models import HeasarcMAXI7YR
from heasarc.models import HeasarcRASS2RXS
from srgcat.models import SkyMaps
from monthplan.models import SurveyHealpixPlate
def get_sources(owner, altsurvey=None):
q = ArtSurveySource.objects.all()
q=q.annotate(exclude=Value(False, output_field=BooleanField()))
try:
params = ArtSurveyParams.objects.get(owner=owner)
except:
return q
if not altsurvey is None:
q=q.filter(survey=altsurvey)
else:
if not params.survey is None:
q=q.filter(survey=params.survey)
if not params.band is None:
q=q.filter(band__slug__exact=params.band)
if not params.glat_min is None:
q=q.filter(bii__gte=params.glat_min)
if not params.glat_max is None:
q=q.filter(bii__lt=params.glat_max)
if not params.ecl_lat_min is None:
q=q.filter(ecl_b__gte=params.ecl_lat_min)
if not params.ecl_lat_max is None:
q=q.filter(ecl_b__lt=params.ecl_lat_max)
if not params.dec_min is None:
q=q.filter(dec__gte=params.dec_min)
if not params.dec_max is None:
q=q.filter(dec__lt=params.dec_max)
if not params.sign_ml_min is None:
""" exclude extended sources from this filtering """
#q=q.filter(Q(ml_sig__gte=params.sign_ml_min) | Q(ext__gte=1))
q=q.filter(ml_sig__gte=params.sign_ml_min)
if not params.sign_ml_max is None:
""" exclude extended sources from this filtering """
#q=q.filter(Q(ml_sig__lt=params.sign_ml_max) | Q(ext__gte=1))
q=q.filter(ml_sig__lt=params.sign_ml_max)
if not params.exposure_min is None:
q=q.filter(exptime__gte=params.exposure_min)
if not params.detlike_min is None:
""" exclude extended sources from this filtering """
q=q.filter(Q(detlike__gte=params.detlike_min) | Q(ext__gte=1))
if not params.detlike_max is None:
""" exclude extended sources from this filtering """
q=q.filter(Q(detlike__lt=params.detlike_max) | Q(ext__gte=1))
if not params.log_nfalse_min is None:
""" exclude extended sources from this filtering """
q=q.filter(Q(log_nfalse__lte=params.log_nfalse_min) | Q(ext__gte=1))
if not params.log_nfalse_max is None:
""" exclude extended sources from this filtering """
q=q.filter(Q(log_nfalse__gt=params.log_nfalse_max) | Q(ext__gte=1))
if not params.log_ml_nfalse_min is None:
""" exclude extended sources from this filtering """
q=q.filter(Q(log_ml_nfalse__lte=params.log_ml_nfalse_min) | Q(ext__gte=1))
if not params.log_ml_nfalse_max is None:
""" exclude extended sources from this filtering """
q=q.filter(Q(log_ml_nfalse__gt=params.log_ml_nfalse_max) | Q(ext__gte=1))
if not params.ext_min is None:
q=q.filter(ext__gte=params.ext_min)
if not params.ext_max is None:
q=q.filter(ext__lt=params.ext_max)
if not params.class_startswith is None:
q=q.filter(object_class__class_name__startswith=params.class_startswith)
if not params.cname_contains is None:
#q=q.filter(cname__startswith=params.cname_startswith)
q=q.filter(cname__contains=params.cname_contains)
if not params.category is None:
q=q.filter(category=params.category)
if not params.exclude_category is None:
q=q.exclude(category=params.exclude_category)
if params.category_unclassified:
q=q.filter(category=None)
if params.gaia_primary:
q=q.filter(gaia3_primary__isnull=False)
if params.turk_possible:
q=q.filter(turkpossible__exact=True)
if params.allwise_primary:
q=q.filter(allwise_primary__isnull=False)
if params.sky == 'allsky':
pass
elif params.sky == 'ru':
q=q.filter(lii__lt=180)
pass
elif params.sky == 'de':
q=q.filter(lii__gt=180)
pass
elif params.sky == 'ru,de':
pass
else:
pass
if params.marshall == None:
pass
elif params.marshall == 'MSFC_FILT':
q=q.filter(ecl_b__gt=82)
pass
elif params.marshall == 'MSFC_EXCL':
q=q.filter(ecl_b__lt=82)
pass
else:
pass
if not params.exclude_survey is None:
#print('params.exclude_survey')
subquery=MetaSource.objects.filter(artsurveysource=OuterRef("pk"))
survey_condition=Q(artsurveysource__survey=params.exclude_survey)
if not params.exclude_band == 'ALL':
#print('params.exclude_band apply E band')
band_condition=Q(artsurveysource__band__slug__exact=params.exclude_band)
else:
""" Take all bands """
#print('params.exclude_band apply ALL bands')
band_condition=Q(artsurveysource__band__slug__contains='E')
print("exclude band: {}".format(params.exclude_band))
""" switched to ML """
if not params.exclude_log_ml_nfalse is None:
#print('params.exclude_log_nfalse apply')
log_nfalse_condition=Q(artsurveysource__log_ml_nfalse__lte=params.exclude_log_ml_nfalse)
else:
""" Take all sources """
#print('params.exclude_log_nfalse default')
log_nfalse_condition=Q(artsurveysource__log_ml_nfalse__lte=-1)
q=q.annotate(exclude=Exists(subquery.filter(survey_condition & band_condition & log_nfalse_condition)))
q=q.exclude(exclude=True)
return q
def get_sources_sorted(owner,refsrc,rmax_deg=10):
ref_crd = SkyCoord(refsrc.ra, refsrc.dec, frame="fk5", unit="deg")
srcs0 = get_sources(owner,altsurvey=refsrc.survey).annotate(separation=Value(0.0, output_field=DecimalField()))
srcs=[]
for src in srcs0:
src_crd = SkyCoord(src.ra, src.dec, frame="fk5", unit="deg")
sep=ref_crd.separation(src_crd).arcmin
if((sep/60) <= rmax_deg):
src.separation=sep
srcs.append(src)
return sorted(srcs, key=lambda x: x.separation, reverse=False)
def get_heasarc_around(owner,refsrc,rmax_arcmin=60):
ref_crd = SkyCoord(refsrc.ra, refsrc.dec, frame="fk5", unit="deg")
srcs0 = get_sources(owner,altsurvey=refsrc.survey).annotate(separation=Value(0.0, output_field=DecimalField()))
srcs=[]
for src in srcs0:
src_crd = SkyCoord(src.ra, src.dec, frame="fk5", unit="deg")
sep=ref_crd.separation(src_crd).arcmin
if(sep <= rmax_arcmin):
heasarcs = src.heasarc.all()
for h in heasarcs:
srcs.append(h)
return srcs
def get_skymap_by_coords(ra, dec):
skymap=None
hp_plate = HEALPix(nside=NSIDE_PLATES,
order=ORDER_PLATES,
frame=FK5())
crd = SkyCoord(ra, dec, frame=FK5(), unit="deg")
healpix_plate = hp_plate.skycoord_to_healpix(crd)
try:
my_hp_plate = SurveyHealpixPlate.objects.get(healpix=healpix_plate)
skymaps = my_hp_plate.skymaps_set.all()
sep_min=10000.0
for smap in skymaps:
smap_crd = SkyCoord(smap.RA_CEN, smap.DE_CEN, frame=FK5(), unit="deg")
sep=crd.separation(smap_crd).arcmin
if(sep < sep_min):
sep_min = sep
skymap = smap
except Exception as e:
print("{}".format(e))
pass
return skymap
def search_by_name_artsurvey(name):
srcs = ArtSurveySource.objects.all().filter(Q(name__contains=name) | Q(cname__contains=name))
return sorted(srcs, key=lambda x: x.name, reverse=False)[:100]
def search_by_coords_artsurvey(ra, dec, rmax=40,maxdist=120.0):
hp = HEALPix(nside=NSIDE_SOURCES, order=ORDER, frame=FK5())
crd = SkyCoord(ra, dec, frame="fk5", unit="deg")
heal = hp.cone_search_skycoord(crd, radius=maxdist*u.arcsecond)
srcs_cone = ArtSurveySource.objects.all().filter(healpix__in=heal)
srcs=[]
for src in srcs_cone:
src_crd = SkyCoord(src.ra, src.dec, frame="fk5", unit="deg")
sep=crd.separation(src_crd).arcsecond
if(sep <= rmax):
src.separation=sep
srcs.append(src)
return sorted(srcs, key=lambda x: x.separation, reverse=False)
#return srcs
def find_heasarc(srcs,maxdist=120.0,minrad=40.0, syserr=0.0, offset=False):
""" version for paper2 """
count=0
hp = HEALPix(nside=NSIDE_SOURCES, order=ORDER, frame=FK5())
excluded_heasarc=0
for s in srcs:
# don't clean at the beginning
#s.heasarc.clear()
crd = SkyCoord(s.ra, s.dec, frame="fk5", unit="deg")
if(offset):
# random angle for estimating spurious matches
alpha = 2 * math.pi * random.random()
position_angle = alpha * 180 / math.pi * u.deg
separation = 45*60 * u.arcsecond
crd=crd.directional_offset_by(position_angle, separation)
heal = hp.cone_search_skycoord(crd, radius=maxdist*u.arcsecond)
try:
heasarcs = HeasarcBase.objects.instance_of(HeasarcSwiftBAT105m,HeasarcXMMSL2,Heasarc4FGL,HeasarcIntegral2020,HeasarcMAXI7YR,HeasarcRASS2RXS)
heasarcs = heasarcs.filter(healpix__in=heal)
#heasarcs = HeasarcBase.objects.all().filter(healpix__in=heal)
#print("{} --> {} HEASARCs".format(s,heasarcs.count()))
except Exception as e:
print('%s' % type(e))
continue
radec_error = sqrt(s.ml_radec_err_98**2 + syserr**2)
heasarcs_clean=heasarcs
if(offset):
for h in heasarcs:
if(h in s.heasarc.all()):
print(s,"*** Already exists in HEASARC list ***",h)
heasarcs_clean=heasarcs_clean.exclude(healpix=h.healpix)
excluded_heasarc+=1
s.heasarc.clear()
if not (heasarcs_clean.count()>0):
continue
#print("{} --> {} HEASARCs === {}".format(s,heasarcs.count(), heasarcs_clean.count()))
for h in heasarcs_clean:
heasarc_crd = SkyCoord(h.ra, h.dec, frame="fk5", unit="deg")
sep=crd.separation(heasarc_crd).arcsecond
if(sep <= radec_error or sep <= (radec_error + h.error_radius)):
s.heasarc.add(h)
s.save()
print("### Total HEASARC excluded {} ###".format(excluded_heasarc))
def find_heasarc_paper1(srcs,maxdist=120.0,minrad=40.0, offset=False):
""" saved version used for paper1 """
count=0
hp = HEALPix(nside=NSIDE_SOURCES, order=ORDER, frame=FK5())
for s in srcs:
s.heasarc.clear()
crd = SkyCoord(s.ra, s.dec, frame="fk5", unit="deg")
if(offset):
# random angle
alpha = 2 * math.pi * random.random()
position_angle = alpha * 180 / math.pi * u.deg
separation = 45*60 * u.arcsecond
crd=crd.directional_offset_by(position_angle, separation)
heal = hp.cone_search_skycoord(crd, radius=maxdist*u.arcsecond)
try:
heasarcs = HeasarcBase.objects.all().filter(healpix__in=heal)
# print("--> HEASARC {}".format(heasarcs.count()))
except:
continue
if(s.ext):
if(s.ext > minrad):
radec_error=float(s.ext)
else:
radec_error=minrad
else:
radec_error=minrad
for h in heasarcs:
heasarc_crd = SkyCoord(h.ra, h.dec, frame="fk5", unit="deg")
sep=crd.separation(heasarc_crd).arcsecond
if(sep <= radec_error or sep <= h.error_radius):
s.heasarc.add(h)
s.save()
def update_marshall(srcs):
count=0
srcs0 = srcs.filter(ecl_b__gt=82.0)
for s in srcs0:
s.marshall=True
s.save()
return srcs0.count()
def update_metasource(srcs,maxdist=120.0,minrad=40.0):
count=0
hp = HEALPix(nside=NSIDE_SOURCES, order=ORDER, frame=FK5())
for s in srcs:
crd = SkyCoord(s.ra, s.dec, frame="fk5", unit="deg")
heal = hp.cone_search_skycoord(crd, radius=maxdist*u.arcsecond)
try:
meta_srcs = MetaSource.objects.all().filter(healpix__in=heal)
except:
continue
sep_min=1000.0
msource=None
for sm in meta_srcs:
meta_crd = SkyCoord(sm.ra, sm.dec, frame="fk5", unit="deg")
sep=crd.separation(meta_crd).arcsecond
if(sep <= sm.radec_error and sep < sep_min):
msource=sm
sep_min=sep
count=count+1
if not (msource):
sep_min=0.0
if(s.ext):
if(s.ext > minrad):
radec_error=float(s.ext)
else:
radec_error=minrad
else:
radec_error=minrad
msource = MetaSource(name=s.name, healpix=s.healpix, ra=s.ra, dec=s.dec, glon=s.lii, glat=s.bii, elon=s.ecl_l, elat=s.ecl_b, radec_error=radec_error)
msource.save()
s.metasource=msource
s.metasource_sep=sep_min
s.save()
return count
def setup_metasource(master_name,slave_name):
try:
srcs1= ArtSurveySource.objects.filter(name__contains=master_name)
print("Master: Found {}".format(srcs1.count()))
for src in srcs1:
print("Master:--> {} --> {}".format(src,src.metasource))
except:
print("Master: Not found: {}".format(master_name))
return
try:
srcs2= ArtSurveySource.objects.filter(name__contains=slave_name)
print("Slave: Found {}".format(srcs2.count()))
for src in srcs2:
print("Slave:--> {} --> {}".format(src,src.metasource))
except:
print("Slave: Not found: {}".format(slave_name))
return
print("*** BEGIN ***")
for master in srcs1:
if not master.metasource:
continue
master_crd = SkyCoord(master.ra, master.dec, frame="fk5", unit="deg")
meta_crd = SkyCoord(master.metasource.ra, master.metasource.dec, frame="fk5", unit="deg")
print("Master {} (<--{})".format(master.metasource,master))
for slave in srcs2:
if not slave.metasource:
continue
if slave.metasource == master.metasource:
print("Already satisfied")
continue
slave_crd = SkyCoord(slave.ra, slave.dec, frame="fk5", unit="deg")
metasource_sep=slave_crd.separation(meta_crd).arcsecond
slave.metasource=master.metasource
slave.metasource_sep=metasource_sep
slave.save()
print("--> Add slave {} at {:.2f} arcseconds".format(slave,metasource_sep))
print("*** END ***")
def get_artsurvey(version, archived=False):
try:
survey = ArtSurvey.objects.get(Q(version=version) & Q(archived=archived))
return survey
except ArtSurvey.DoesNotExist:
return None
def get_energy_band(slug):
try:
band = EnergyBand.objects.get(slug=slug)
return band
except EnergyBand.DoesNotExist:
return None
def find_skymaps():
skymaps = SkyMaps.objects.all()
print("Selected {} skymaps".format(skymaps.count()))
srcs = ArtSurveySource.objects.all().filter(SMAPNR__exact=0)
print("Selected {} total".format(srcs.count()))
for src in srcs:
crd = SkyCoord(src.ra, src.dec, frame="fk5", unit="deg")
min_sep=400
for smap in skymaps:
crd_map = SkyCoord(smap.RA_CEN, smap.DE_CEN, frame="fk5", unit="deg")
sep=crd.separation(crd_map).deg
if(sep < min_sep):
min_sep=sep
src.SMAPNR=smap.SMAPNR
if(min_sep < 360):
print("Set SkyMap for {} at offset {} degrees".format(src,min_sep))
src.save()
def select_turkpossible(survey_version,slug):
minrad=40
maxdist=600
try:
band = EnergyBand.objects.get(slug=slug)
except EnergyBand.DoesNotExist:
print("EnergyBand not found")
return
"""
try:
survey = ArtSurvey.objects.get(version=survey_version)
print("use ArtSurvey {}".format(survey))
except ArtSurvey.DoesNotExist:
print("ArtSurvey {} not found".format(survey_version))
return
srcs = survey.artsurveysource_set.all().filter(band=band)
"""
srcs = ArtSurveySource.objects.all().filter(Q(survey__version=survey_version)&Q(band=band))
print("Selected {}".format(srcs.count()))
for src in srcs:
src.turkpossible=True
src.save()
for src in srcs:
heasarcs = src.heasarc.all()
if (heasarcs.filter(instance_of=HeasarcSwiftBAT105m).count()):
src.turkpossible=False
src.save()
if (heasarcs.filter(instance_of=HeasarcIntRefCat).count()):
src.turkpossible=False
src.save()
if (heasarcs.filter(instance_of=HeasarcXTEASSCAT).count()):
src.turkpossible=False
src.save()
selection=srcs.filter(turkpossible__exact=True)
print("Total selected {}".format(selection.count()))
def ps1_getimages(ra,dec,size=240,filters="grizy"):
"""Query ps1filenames.py service to get a list of images
ra, dec = position in degrees
size = image size in pixels (0.25 arcsec/pixel)
filters = string with filters to include
Returns a table with the results
"""
service = "https://ps1images.stsci.edu/cgi-bin/ps1filenames.py"
url = ("{service}?ra={ra}&dec={dec}&size={size}&format=fits"
"&filters={filters}").format(**locals())
table = Table.read(url, format='ascii')
return table
def ps1_geturl(ra, dec, size=240, output_size=None, filters="grizy", format="jpg", color=False):
"""Get URL for images in the table
ra, dec = position in degrees
size = extracted image size in pixels (0.25 arcsec/pixel)
output_size = output (display) image size in pixels (default = size).
output_size has no effect for fits format images.
filters = string with filters to include
format = data format (options are "jpg", "png" or "fits")
color = if True, creates a color image (only for jpg or png format).
Default is return a list of URLs for single-filter grayscale images.
Returns a string with the URL
"""
if color and format == "fits":
raise ValueError("color images are available only for jpg or png formats")
if format not in ("jpg","png","fits"):
raise ValueError("format must be one of jpg, png, fits")
table = ps1_getimages(ra,dec,size=size,filters=filters)
url = ("https://ps1images.stsci.edu/cgi-bin/fitscut.cgi?"
"ra={ra}&dec={dec}&size={size}&format={format}").format(**locals())
if output_size:
url = url + "&output_size={}".format(output_size)
# sort filters from red to blue
flist = ["yzirg".find(x) for x in table['filter']]
table = table[np.argsort(flist)]
if color:
if len(table) > 3:
# pick 3 filters
table = table[[0,len(table)//2,len(table)-1]]
for i, param in enumerate(["red","green","blue"]):
url = url + "&{}={}".format(param,table['filename'][i])
else:
urlbase = url + "&red="
url = []
for filename in table['filename']:
url.append(urlbase+filename)
return url
def skyview_query(ra, dec, size=240, survey=None):
crd = SkyCoord(ra, dec, frame="fk5", unit="deg")
urls = SkyView.get_image_list(position=crd, survey=[survey])
r = requests.get(urls[0])
return r.content
def ps1_getcolorim(ra, dec, size=240, output_size=None, filters="grizy", format="jpg"):
"""Get color image at a sky position
ra, dec = position in degrees
size = extracted image size in pixels (0.25 arcsec/pixel)
output_size = output (display) image size in pixels (default = size).
output_size has no effect for fits format images.
filters = string with filters to include
format = data format (options are "jpg", "png")
Returns the image
"""
if format not in ("jpg","png"):
raise ValueError("format must be jpg or png")
url = ps1_geturl(ra,dec,size=size,filters=filters,output_size=output_size,format=format,color=True)
r = requests.get(url)
#im = Image.open(BytesIO(r.content))
return r.content
def ps1_getgrayim(ra, dec, size=240, output_size=None, filter="g", format="jpg"):
"""Get grayscale image at a sky position
ra, dec = position in degrees
size = extracted image size in pixels (0.25 arcsec/pixel)
output_size = output (display) image size in pixels (default = size).
output_size has no effect for fits format images.
filter = string with filter to extract (one of grizy)
format = data format (options are "jpg", "png")
Returns the image
"""
if format not in ("jpg","png"):
raise ValueError("format must be jpg or png")
if filter not in list("grizy"):
raise ValueError("filter must be one of grizy")
url = ps1_geturl(ra,dec,size=size,filters=filter,output_size=output_size,format=format)
r = requests.get(url[0])
#im = Image.open(BytesIO(r.content))
return r.content
def ps1_fitsim(ra, dec, size=240, output_size=None, filter="g", format="fits"):
if format not in ("fits"):
raise ValueError("format must be fits")
if filter not in list("grizy"):
raise ValueError("filter must be one of grizy")
url = ps1_geturl(ra,dec,size=size,filters=filter,output_size=output_size,format=format)
r = requests.get(url[0])
#im = Image.open(BytesIO(r.content))
return r.content