Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MOC calculation #154

Merged
merged 4 commits into from
May 28, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 2 additions & 10 deletions gwemopt/args.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,6 @@ def parse_args(args):
parser.add_argument("--doPerturbativeTiling", action="store_true", default=False)
parser.add_argument("--doOrderByObservability", action="store_true", default=False)

parser.add_argument("--doUseCatalog", action="store_true", default=False)
parser.add_argument("--catalogDir", help="catalog directory", default=CATALOG_DIR)
parser.add_argument(
"--catalog", help="Galaxy catalog name (e.g GLADE)", default=None
Expand All @@ -54,7 +53,6 @@ def parse_args(args):
type=int,
default=2000,
)
parser.add_argument("--doObservability", action="store_true", default=False)

parser.add_argument("--doEfficiency", action="store_true", default=False)
parser.add_argument(
Expand Down Expand Up @@ -83,7 +81,6 @@ def parse_args(args):

parser.add_argument("--Ndet", default=1, type=int)
parser.add_argument("--nside", default=256, type=int)
parser.add_argument("--DScale", default=1.0, type=float)
parser.add_argument("--Tobs", default="0.0,1.0")

parser.add_argument("--mindiff", default=0.0, type=float)
Expand Down Expand Up @@ -111,8 +108,6 @@ def parse_args(args):
parser.add_argument("--doSingleExposure", action="store_true", default=False)
parser.add_argument("--filters", default="r,g,r")
parser.add_argument("--doAlternatingFilters", action="store_true", default=False)
parser.add_argument("--doRASlices", action="store_true", default=False)
parser.add_argument("--nside_down", default=2, type=int)
parser.add_argument("--max_filter_sets", default=4, type=int)
parser.add_argument("--iterativeOverlap", default=0.0, type=float)
parser.add_argument("--maximumOverlap", default=1.0, type=float)
Expand All @@ -137,10 +132,8 @@ def parse_args(args):
parser.add_argument("--doBlocks", action="store_true", default=False)
parser.add_argument("--Nblocks", default=4, type=int)

parser.add_argument("--doRASlice", action="store_true", default=False)
parser.add_argument("--raslice", default="0.0,24.0")

parser.add_argument("--absmag", default=-15.0, type=float)
parser.add_argument("--raslice", default=None)

parser.add_argument(
"--galactic_limit", help="Galactic limit.", default=0.0, type=float
Expand All @@ -155,10 +148,9 @@ def parse_args(args):
parser.add_argument("--true_inclination", default=0.0, type=float)

parser.add_argument("--doTrueLocation", action="store_true", default=False)
parser.add_argument("--observability_thresh", default=0.05, type=float)
parser.add_argument("--doObservabilityExit", action="store_true", default=False)

parser.add_argument("--inclination", action="store_true", default=False)
parser.add_argument("--projection", default="astro mollweide")

parser.add_argument("--solverType", default="heuristic")
parser.add_argument("--milpSolver", default="PULP_CBC_CMD")
Expand Down
45 changes: 28 additions & 17 deletions gwemopt/catalogs/get.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,10 @@

import healpy as hp
import numpy as np
from astropy.table import Table
from astroquery.vizier import Vizier
from ligo.skymap import distance
from ligo.skymap.bayestar import derasterize
from scipy.stats import norm

from gwemopt.catalogs.clu import CluCatalog
Expand Down Expand Up @@ -63,7 +65,7 @@ def get_catalog(params, map_struct, export_catalog: bool = True):

n, cl = params["powerlaw_n"], params["powerlaw_cl"]

prob_scaled = copy.deepcopy(map_struct["prob"])
prob_scaled = copy.deepcopy(map_struct["skymap_raster"]["PROB"])
prob_sorted = np.sort(prob_scaled)[::-1]
prob_indexes = np.argsort(prob_scaled)[::-1]
prob_cumsum = np.cumsum(prob_sorted)
Expand All @@ -72,34 +74,42 @@ def get_catalog(params, map_struct, export_catalog: bool = True):
prob_scaled = prob_scaled**n

ipix = hp.ang2pix(
map_struct["nside"],
params["nside"],
np.array(cat_df["ra"]),
np.array(cat_df["dec"]),
lonlat=True,
)

if "distnorm" in map_struct:
if map_struct["distnorm"] is not None:
if "DISTNORM" in map_struct:
if map_struct["skymap_raster"]["DISTNORM"] is not None:
# creat an mask to cut at 3 sigma in distance
mask = np.zeros(len(cat_df["distmpc"]))

# calculate the moments from distmu, distsigma and distnorm
mom_mean, mom_std, mom_norm = distance.parameters_to_moments(
map_struct["distmu"], map_struct["distsigma"]
)

condition_indexer = np.where(
(cat_df["distmpc"] < (mom_mean[ipix] + (3 * mom_std[ipix])))
& (cat_df["distmpc"] > (mom_mean[ipix] - (3 * mom_std[ipix])))
(
cat_df["distmpc"]
< (
map_struct["skymap_raster"]["DISTMEAN"][ipix]
+ (3 * map_struct["skymap_raster"]["DISTSTD"][ipix])
)
)
& (
cat_df["distmpc"]
> (
map_struct["skymap_raster"]["DISTMEAN"][ipix]
- (3 * map_struct["skymap_raster"]["DISTSTD"][ipix])
)
)
)
mask[condition_indexer] = 1

s_loc = (
prob_scaled[ipix]
* (
map_struct["distnorm"][ipix]
map_struct["skymap_raster"]["DISTNORM"][ipix]
* norm(
map_struct["distmu"][ipix], map_struct["distsigma"][ipix]
map_struct["skymap_raster"]["DISTMU"][ipix],
map_struct["skymap_raster"]["DISTSIGMA"][ipix],
).pdf(cat_df["distmpc"])
)
** params["powerlaw_dist_exp"]
Expand Down Expand Up @@ -202,7 +212,7 @@ def get_catalog(params, map_struct, export_catalog: bool = True):
s_mass = s_loc * (1 + (alpha_mass * beta_mass * s_mass))

s = np.array(s_loc * Slum * sdet)
prob = np.zeros(map_struct["prob"].shape)
prob = np.zeros(map_struct["skymap_raster"]["PROB"].shape)
if params["galaxy_grade"] == "Sloc":
for j in range(len(ipix)):
prob[ipix[j]] += s_loc[j]
Expand All @@ -220,11 +230,12 @@ def get_catalog(params, map_struct, export_catalog: bool = True):
"You are trying to use a galaxy grade that is not implemented yet."
)

prob[np.isnan(prob)] = 0.0
prob = prob / np.sum(prob)

map_struct["prob_catalog"] = prob
if params["doUseCatalog"]:
map_struct["prob"] = prob
map_struct["skymap_raster"]["PROB"] = prob
skymap = map_struct["skymap_raster"].copy()
map_struct["skymap"] = derasterize(skymap)

cat_df["grade"] = grade
cat_df["S"] = s
Expand Down
4 changes: 2 additions & 2 deletions gwemopt/chipgaps/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,5 +2,5 @@
Module for chip gaps.
"""

from gwemopt.chipgaps.decam import get_decam_quadrant_ipix
from gwemopt.chipgaps.ztf import get_ztf_quadrant_ipix
from gwemopt.chipgaps.decam import get_decam_quadrant_moc
from gwemopt.chipgaps.ztf import get_ztf_quadrant_moc
17 changes: 10 additions & 7 deletions gwemopt/chipgaps/decam.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
from astropy import units as u
from astropy.coordinates import ICRS, SkyCoord
from astropy.wcs import WCS
from mocpy import MOC


class DECamtile:
Expand Down Expand Up @@ -210,7 +211,7 @@ def get_decam_ccds(ra, dec, save_footprint=False):
return np.transpose(offsets, (2, 0, 1))


def get_decam_quadrant_ipix(nside, ra, dec):
def get_decam_quadrant_moc(ra, dec):
ccd_coords = get_decam_ccds(0, 0)

skyoffset_frames = SkyCoord(ra, dec, unit=u.deg).skyoffset_frame()
Expand All @@ -219,13 +220,15 @@ def get_decam_quadrant_ipix(nside, ra, dec):
unit=u.deg,
frame=skyoffset_frames[:, np.newaxis, np.newaxis],
).transform_to(ICRS)
ccd_xyz = np.moveaxis(ccd_coords_icrs.cartesian.xyz.value, 0, -1)[0]

ipixs = []
for xyz in ccd_xyz:
ipix = hp.query_polygon(nside, xyz)
ipixs.append(ipix.tolist())
return ipixs
moc = None
for subfield_id, coords in enumerate(ccd_coords_icrs[0]):
if moc is None:
moc = MOC.from_polygon_skycoord(coords)
else:
moc = moc + MOC.from_polygon_skycoord(coords)

return moc


def ccd_xy_to_radec(alpha_p, delta_p, ccd_centers_xy):
Expand Down
18 changes: 11 additions & 7 deletions gwemopt/chipgaps/ztf.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
from astropy import units as u
from astropy.coordinates import ICRS, SkyCoord
from astropy.wcs import WCS
from mocpy import MOC


class ZTFtile:
Expand Down Expand Up @@ -273,7 +274,7 @@ def get_ztf_quadrants():
return np.transpose(offsets, (2, 0, 1))


def get_ztf_quadrant_ipix(nside, ra, dec, subfield_ids=None):
def get_ztf_quadrant_moc(ra, dec, subfield_ids=None):
quadrant_coords = get_ztf_quadrants()

skyoffset_frames = SkyCoord(ra, dec, unit=u.deg).skyoffset_frame()
Expand All @@ -282,13 +283,16 @@ def get_ztf_quadrant_ipix(nside, ra, dec, subfield_ids=None):
unit=u.deg,
frame=skyoffset_frames[:, np.newaxis, np.newaxis],
).transform_to(ICRS)
quadrant_xyz = np.moveaxis(quadrant_coords_icrs.cartesian.xyz.value, 0, -1)[0]

ipixs = []
for subfield_id, xyz in enumerate(quadrant_xyz):
moc = None
for subfield_id, coords in enumerate(quadrant_coords_icrs[0]):
if subfield_ids is not None:
if subfield_id not in subfield_ids:
continue
ipix = hp.query_polygon(nside, xyz)
ipixs.append(ipix.tolist())
return ipixs

if moc is None:
moc = MOC.from_polygon_skycoord(coords)
else:
moc = moc + MOC.from_polygon_skycoord(coords)

return moc
Loading
Loading