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.

Sign up for GitHub

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

io.skymap.read_skymap(): no_hdu parameter #159

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

io.skymap.read_skymap(): no_hdu parameter #159

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
be17053
add no_hdu parameter to read_skymap method, to avoid creating the rat…
Theodlz Dec 2, 2024
3d53431
single quotes to double quotes for consistency
Theodlz Dec 2, 2024
2d74845
remove no_hdu parameter, and replace by no_plotting, no_catalog, and …
Theodlz Dec 3, 2024
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
149 changes: 87 additions & 62 deletions gwemopt/io/skymap.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -214,94 +214,105 @@ def read_inclination(skymap, params, map_struct):
return map_struct


def read_skymap(params, map_struct=None):
def read_skymap(params, map_struct=None, no_plotting=False, no_efficiency=False, no_catalog=False):
"""
Read in a skymap and return a map_struct

:param params: dictionary of parameters
:param map_struct: dictionary of map parameters
:return: map_struct
:param no_hdu: do not create an HDU, default False
:return: params, map_struct
"""

geometry = params["geometry"]
if geometry is not None:
if geometry == "2d":
params["do_3d"] = False
else:
params["do_3d"] = True
# some of the variables that this function creates
# are only necessary if we plan on doing:
# 1. loading catalogs (no_catalog=False)
# 2. plotting (no_plotting=False)
# 3. computing efficiency (no_efficiency=False)
# if we don't plan on doing some of these, we can skip
# some of the calculations to save time and memory

if params.get("geometry") is None:
pass
elif params.get("geometry") == "2d":
params["do_3d"] = False
else:
params["do_3d"] = True

if map_struct is None and params.get("skymap") is None:
raise ValueError("No skymap provided")

if map_struct is None:
# Let's just figure out what's in the skymap first
skymap_path = params["skymap"]

params["name"] = Path(skymap_path).stem

is_3d = False
t_obs = Time.now()
do_3d = False

with fits.open(skymap_path) as hdul:
for x in hdul:
if "DATE-OBS" in x.header:
t_obs = Time(x.header["DATE-OBS"], format="isot")

elif "EVENTMJD" in x.header:
t_obs_mjd = x.header["EVENTMJD"]
t_obs = Time(t_obs_mjd, format="mjd")
t_obs = Time(x.header["EVENTMJD"], format="mjd")

if ("DISTMEAN" in x.header) | ("DISTSTD" in x.header):
is_3d = True
do_3d = True

# Set GPS time from skymap, if not specified. Defaults to today
params["eventtime"] = t_obs
if params["gpstime"] is None:
params["gpstime"] = t_obs.gps

if "do_3d" not in params:
params["do_3d"] = is_3d
params["do_3d"] = do_3d

map_struct = {}

filename = params["skymap"]

if params["do_3d"]:
skymap = read_sky_map(filename, moc=True, distances=True)

if "PROBDENSITY_SAMPLES" in skymap.columns:
if params["inclination"]:
map_struct = read_inclination(skymap, params, map_struct)

skymap.remove_columns(
[
f"{name}_SAMPLES"
for name in [
"PROBDENSITY",
"DISTMU",
"DISTSIGMA",
"DISTNORM",
]
skymap = read_sky_map(skymap_path, moc=True, distances=params["do_3d"])

if (
params["do_3d"]
and "PROBDENSITY_SAMPLES" in skymap.columns
):
if params["inclination"]:
map_struct = read_inclination(skymap, params, map_struct)

skymap.remove_columns(
[
f"{name}_SAMPLES"
for name in [
"PROBDENSITY",
"DISTMU",
"DISTSIGMA",
"DISTNORM",
]
)
]
)

map_struct["skymap"] = skymap
else:
skymap = read_sky_map(filename, moc=True, distances=False)
map_struct["skymap"] = skymap
map_struct["skymap"] = skymap

level, ipix = ah.uniq_to_level_ipix(map_struct["skymap"]["UNIQ"])
nside = ah.level_to_nside(level)
ra, dec = ah.healpix_to_lonlat(ipix, nside, order="nested")
map_struct["skymap"]["ra"] = ra.deg
map_struct["skymap"]["dec"] = dec.deg

map_struct["skymap_raster"] = rasterize(
map_struct["skymap"], order=hp.nside2order(int(params["nside"]))
)
if not (no_efficiency and no_plotting):
# rasterize the skymap, if plotting or efficiency calculations are needed
map_struct["skymap_raster"] = rasterize(
map_struct["skymap"], order=hp.nside2order(int(params["nside"]))
)

peak = map_struct["skymap_raster"][
map_struct["skymap_raster"]["PROB"]
== np.max(map_struct["skymap_raster"]["PROB"])
]
map_struct["center"] = SkyCoord(peak["ra"][0] * u.deg, peak["dec"][0] * u.deg)
if not no_plotting:
# find the peak of the skymap, if plotting is needed
peak = map_struct["skymap_raster"][
map_struct["skymap_raster"]["PROB"]
== np.max(map_struct["skymap_raster"]["PROB"])
]
map_struct["center"] = SkyCoord(peak["ra"][0] * u.deg, peak["dec"][0] * u.deg)

map_struct["skymap_schedule"] = map_struct["skymap"].copy()

Expand All @@ -321,9 +332,16 @@ def read_skymap(params, map_struct=None):
cumprob = np.cumsum(prob)
ii = np.where(cumprob > params["confidence_level"])[0]
map_struct["skymap_raster_schedule"]["PROB"][ind[ii]] = 0.0
map_struct["skymap_schedule"] = derasterize(map_struct["skymap_raster"].copy())
map_struct["skymap_schedule"] = derasterize(map_struct["skymap_raster_schedule"].copy())

if "DISTMU" in map_struct["skymap_raster"].columns:
if no_catalog:
# remove rasterized skymap schedule, if loading a catalog is not needed
del map_struct["skymap_raster_schedule"]

if (
"skymap_raster" in map_struct
and "DISTMU" in map_struct["skymap_raster"].columns
):
(
map_struct["skymap_raster"]["DISTMEAN"],
map_struct["skymap_raster"]["DISTSTD"],
Expand All @@ -333,20 +351,27 @@ def read_skymap(params, map_struct=None):
map_struct["skymap_raster"]["DISTSIGMA"],
)

extra_header = [
("PIXTYPE", "HEALPIX", "HEALPIX pixelisation"),
("ORDERING", "NESTED", "Pixel ordering scheme: RING, NESTED, or NUNIQ"),
("COORDSYS", "C", "Ecliptic, Galactic or Celestial (equatorial)"),
(
"MOCORDER",
moc.uniq2order(map_struct["skymap"]["UNIQ"].max()),
"MOC resolution (best order)",
),
("INDXSCHM", "EXPLICIT", "Indexing: IMPLICIT or EXPLICIT"),
]

hdu = fits.table_to_hdu(map_struct["skymap_raster"])
hdu.header.extend(extra_header)
map_struct["hdu"] = hdu
map_struct["hdu"] = None
if not no_plotting:
# create an HDU for the skymap, if plotting is needed
extra_header = [
("PIXTYPE", "HEALPIX", "HEALPIX pixelisation"),
("ORDERING", "NESTED", "Pixel ordering scheme: RING, NESTED, or NUNIQ"),
("COORDSYS", "C", "Ecliptic, Galactic or Celestial (equatorial)"),
(
"MOCORDER",
moc.uniq2order(map_struct["skymap"]["UNIQ"].max()),
"MOC resolution (best order)",
),
("INDXSCHM", "EXPLICIT", "Indexing: IMPLICIT or EXPLICIT"),
]

hdu = fits.table_to_hdu(map_struct["skymap_raster"])
hdu.header.extend(extra_header)
map_struct["hdu"] = hdu

if no_efficiency:
# remove the rasterized skymap, if efficiency calculations are not needed
del map_struct["skymap_raster"]

return params, map_struct
Loading