Select Git revision
MaskCreator.py 12.84 KiB
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Feb 3 14:03:21 2020
@author: mmayer
"""
from astropy.io import fits #import packages
import numpy as np
import astropy.units as u
from astropy.coordinates import Angle
from astropy.wcs import WCS as PyWCS
def GeneralMaskCreator(Directory, skyfield='', config='', version='', prefix='',
srcregion='',backregion='',exclregion='',
catfile='',phirange='',
templatefile='',templatefilteredfile='',
Rate0 = 3E-1, Radius0 = 0.5, MinLike = 50, extended=False,
):
"""
This method maps ds9 regions to fits mask files, optionally excluding sources in a catalog.
The following parameters are important:
- Directory: The directory in which all files are
- srcregion, backregion: The names (without ".reg") of the source and background regions in the directory.
- phirange: For circles & ellipse: Which ranges of angles to include (--> To make "pizza shaped regions")
- catfile: Path to the catalog which should be excluded (e.g. eRASS1 catalog)
- templatefile: Template image, whose coordinate system and pixel size will be used to create fits masks
- Rate0, Radius0: Scaling of exclusion radius as Radius = Radius0 * sqrt(Rate/Rate0), Radius0 in units of arcmin
- MinLike: Sources above this likelihood will be masked
- extended: Whether to also mask extended sources
"""
if catfile:
Cat = fits.open(catfile)
DATA = Cat[1].data
# Names = DATA.names
DETLIKE = DATA['DET_LIKE_0']
EXTLIKE = DATA['EXT_LIKE']
if extended:
EroMask = (DETLIKE >= MinLike)
else:
EroMask = (DETLIKE >= MinLike) * (EXTLIKE == 0)
FILT_DATA = DATA[EroMask]
DEC = FILT_DATA.field('DEC')
RA = FILT_DATA.field('RA')
RATE = FILT_DATA.field('ML_RATE_0')
CTS = FILT_DATA.field('ML_CTS_0')
# Radii = Radius0 * (CTS/Counts0)**0.5 * 60 ### TO ARCSEC
Radii = Radius0 * (RATE/Rate0)**0.5 * 60 ### TO ARCSEC
if extended:
EXTENT = FILT_DATA['EXT'] ### in arcsec!
ExtRadii = EXTENT*2
Radii = np.max([Radii,ExtRadii], axis=0)
""" FIRST, OPEN TEMPLATE IMAGE FILE, TO GET STRUCTURE CORRECTLY! """
if templatefile:
file = templatefile
else:
file = '%s/%s_%s_024_EventList_%s_c%s.fits' % (Directory, prefix, skyfield, config, version) ### TEMPLATEFILE! ####
# with fits.open(file) as hdul:
hdul = fits.open(file)
newhdul = hdul[:1]
header = newhdul[0].header
PixSize = abs(3600*header['CDELT1']) #IN ARCSEC
Image = newhdul[0].data
# Mask = np.zeros_like(Image)
size=Image.shape
w = PyWCS(header)
xmin = 0
xmax = size[1]
ymin = 0
ymax = size[0]
xx = np.arange(xmin,xmax)
yy = np.arange(ymin,ymax)
X, Y = np.meshgrid(xx,yy)
# Y, X = np.meshgrid(yy,xx)
# print(X.shape,Image.shape)
CatMask = np.ones_like(X)
if catfile:
ramax1, decmin1 = w.all_pix2world(xmin,ymin,0)
ramin1, decmax1 = w.all_pix2world(xmax,ymax,0)
ramax2, decmax2 = w.all_pix2world(xmin,ymax,0)
ramin2, decmin2 = w.all_pix2world(xmax,ymin,0)
ramin, decmin = min(ramin1, ramin2), min(decmin1, decmin2)
ramax, decmax = max(ramax1, ramax2), max(decmax1, decmax2)
eroMask = (RA >= ramin) * (RA <= ramax) * (DEC >= decmin) * (DEC <= decmax)
# Rero = ((ra-RA)**2 * np.cos(dec/180*np.pi)**2 + (dec-DEC)**2)**0.5 * 3600
# eroMask = (Rero < rmax)
LocRA, LocDEC, LocRadii = (RA[eroMask], DEC[eroMask], Radii[eroMask])
xero,yero = w.all_world2pix(LocRA,LocDEC,0)
rpix = np.max([1*np.ones_like(LocRadii), LocRadii/PixSize],axis=0)
# print('Pixel radii',rpix)
for j in range(len(LocRA)):
D = ((xero[j]-X)**2+(yero[j]-Y)**2)**0.5
locmask = (D <= rpix[j])
CatMask *= (1 - locmask)
""" OPEN DS9 REGION FILE TO SEE WHAT WERE DEALING WITH """
if srcregion == '':
if exclregion:
Regions = [exclregion]
else:
Regions = []
elif backregion != '': ### DO WE WANT A BACKGROUND REGION TOO?
Regions = [srcregion,backregion]
else:
Regions = [srcregion]
for i in range(len(Regions)):
region = Regions[i]
regfile = open('%s/%s.reg' % (Directory, region),'r')
start = 0
Lines = []
TotMask = np.zeros_like(X)
for line in regfile:
if start == 1:
# Line = line
# break
Lines.append(line)
if 'fk5' in line:
start = 1
for Line in Lines: ### Allows multiple shapes to be combined in a single region file!
### HAVE THE LINE NOW! ###
mode,content = Line.split('(')
if mode == 'ellipse': # Figure out if ellipse or annulus
if len(content.split(')')[0].split(',')) == 5:
pass
if len(content.split(')')[0].split(',')) == 7:
mode = 'ellannulus'
### Convert coordinates to image and create the mask
if mode in ('circle','ellipse'):
if mode == 'ellipse':
ra,dec,a,b,phi0 = content.split(')')[0].split(',')
if mode == 'circle':
ra,dec,r = content.split(')')[0].split(',')
a, b = r, r
phi0 = 180
ra = Angle(ra+' hours').deg
dec = Angle(dec,unit=u.deg).deg
a = Angle(a).arcsec
b = Angle(b).arcsec
phi0 = float(phi0)
rmax = max(a,b) ### Radius at which to start excluding Sources
center = w.all_world2pix(ra,dec,0)
a0,b0 = (a/PixSize, b/PixSize)### semiaxes
D = ((center[0]-X)**2+(center[1]-Y)**2)**0.5
Dir = np.arctan2((Y-center[1]),(X-center[0]))*180/np.pi
DAngle = (Dir-phi0)/180*np.pi
r0 = a0*b0/(b0**2*np.cos(DAngle)**2+a0**2*np.sin(DAngle)**2)**0.5
Mask = (D <= r0) ## Basic Mask
if phirange: ### LIMIT range, to create complex sectors of elliptical annuli!
phimin, phimax = phirange
if phimax > phimin: ### NORMAL ORDERING!
ModDir = Dir % 360
AngleMask = (ModDir >= phimin) * (ModDir <= phimax)
else: ### CROSSING THE ZEROPOINT!
ModDir = Dir % 360
AngleMask = ((ModDir >= phimin) + (ModDir <= phimax))
Mask = Mask * AngleMask
elif mode == 'ellannulus':### Slightly different for elliptical annulus
ra,dec,a0,b0,a1,b1,phi0 = content.split(')')[0].split(',')
ra = Angle(ra+' hours').deg
dec = Angle(dec,unit=u.deg).deg
a0 = Angle(a0).arcsec
b0 = Angle(b0).arcsec
a1 = Angle(a1).arcsec
b1 = Angle(b1).arcsec
phi0 = float(phi0)
rmax = max(a1,b1) ### Radius at which to start excluding Sources
center = w.all_world2pix(ra,dec,0)
a0,b0 = (a0/PixSize, b0/PixSize)### Inner semiaxes
a1,b1 = (a1/PixSize, b1/PixSize)### Outer semiaxes
D = ((center[0]-X)**2+(center[1]-Y)**2)**0.5
Dir = np.arctan2((Y-center[1]),(X-center[0]))*180/np.pi
DAngle = (Dir-phi0)/180*np.pi
r0 = a0*b0/(b0**2*np.cos(DAngle)**2+a0**2*np.sin(DAngle)**2)**0.5
r1 = a1*b1/(b1**2*np.cos(DAngle)**2+a1**2*np.sin(DAngle)**2)**0.5
Mask = (D >= r0)*(D <= r1) ## Basic Mask
if phirange: ### LIMIT range, to create complex sectors of elliptical annuli!
phimin, phimax = phirange
if phimax > phimin: ### NORMAL ORDERING!
ModDir = Dir % 360
AngleMask = (ModDir >= phimin) * (ModDir <= phimax)
else: ### CROSSING THE ZEROPOINT!
ModDir = Dir % 360
AngleMask = ((ModDir >= phimin) + (ModDir <= phimax))
Mask = Mask * AngleMask
if mode == 'box':
ra,dec,a,b,phi0 = content.split(')')[0].split(',')
ra = Angle(ra+' hours').deg
dec = Angle(dec,unit=u.deg).deg
a = Angle(a).arcsec
b = Angle(b).arcsec
phi0 = float(phi0)
rmax = max(a,b) ### Radius at which to start excluding Sources
center = w.all_world2pix(ra,dec,0)
a0,b0 = (a/PixSize, b/PixSize)### semiaxes
RelX, RelY = X - center[0], Y - center[1]
ROT = np.array([[np.cos(phi0/180*np.pi), np.sin(phi0/180*np.pi)], ### NEGATIE ROTATION
[-np.sin(phi0/180*np.pi), np.cos(phi0/180*np.pi)]])
# REL = np.transpose([RelX,RelY],axes=(1,2,0))
# M = np.array([[np.dot(ROT,j) for j in i] for i in REL])
# RotX, RotY = np.transpose(M,axes=(2,0,1))
REL = np.transpose([RelX,RelY],axes=(1,0,2))
M = np.dot(ROT,REL)
RotX, RotY = M
Mask = (abs(RotX) <= a0/2) * (abs(RotY) <= b0/2)
elif mode == 'polygon':
from shapely.geometry import Point, Polygon
radeclist = [float(j) for j in content.split(')')[0].split(',')]
ra, dec = np.reshape(radeclist, (-1,2)).T
x, y = w.all_world2pix(ra,dec,0)
minx, maxx, miny, maxy = int(np.min(x)), int(np.max(x))+1, int(np.min(y)), int(np.max(y))+1
Points = np.transpose(np.meshgrid(np.arange(minx,maxx), np.arange(miny,maxy)), (1,2,0))
PointCoords = np.concatenate(Points)
Points = [Point(p) for p in PointCoords]
Poly = Polygon(np.transpose([x,y]))
mask = np.array([p.within(Poly) for p in Points])
MaskPoints = PointCoords[mask]
Mask = np.zeros_like(CatMask)
for p in MaskPoints:
Mask[p[1], p[0]] = 1
### NOW INCLUDE SOURCE CATLOG! ###
if catfile and srcregion:
Mask = Mask * CatMask
TotMask += Mask ### ITERATIVELY!
if not srcregion and exclregion:
InclMask = (TotMask == 0)
continue
TotMask = TotMask.astype(int)
# print('Old',np.min(newhdul.data),np.max(newhdul.data))
# print('Mask',np.min(TotMask),np.max(TotMask))
newhdul[0].data = TotMask
# print('New',np.min(newhdul[0].data),np.max(newhdul[0].data))
NewName = '%s/%s_mask.fits' % (Directory, region)
newhdul.writeto(NewName,overwrite=True)
if len(srcregion):
return
TotMask = CatMask.astype(int)
if exclregion:
TotMask *= InclMask
MaskedImage = Image*TotMask
newhdul[0].data = TotMask.astype(np.int32)
NewName = '%s/%s_mask.fits' % (Directory, 'pure')
newhdul.writeto(NewName,overwrite=True)
newhdul[0].data = MaskedImage
if templatefilteredfile:
NewName = templatefilteredfile
else:
NewName = '%s/%s_filtered.fits' % (Directory, 'pure')
newhdul.writeto(NewName,overwrite=True)
return
""" Example call: """
GeneralMaskCreator('Vela', srcregion='TestPolygon', backregion='Back', #directory and file paths of region file
templatefile='Vela/sm04_merged_020_WideImage_c020.fits', #Image encompassing region to be masked
catfile='CrossMatching/final_e1_SourceCat1B_211229_poscorr_mpe_clean.fits', #path to catalog file
Rate0 = 3E-1, Radius0 = 0.5, #exclusion radius of 30" at 0.3 ct/s
MinLike = 10,extended=False, #minimum likelihood of DETLIKE_0 = 10, no extended sources are masked
)