From be17053039586d882563c34fb310a71266bf1e45 Mon Sep 17 00:00:00 2001 From: Theodlz Date: Mon, 2 Dec 2024 15:17:16 -0800 Subject: [PATCH 1/3] add no_hdu parameter to read_skymap method, to avoid creating the rather heavy hdu variable if not necessarily (plotting only) --- gwemopt/io/skymap.py | 109 ++++++++++++++++++++++--------------------- 1 file changed, 55 insertions(+), 54 deletions(-) diff --git a/gwemopt/io/skymap.py b/gwemopt/io/skymap.py index bb39444..8ae31b9 100644 --- a/gwemopt/io/skymap.py +++ b/gwemopt/io/skymap.py @@ -214,21 +214,25 @@ 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_hdu=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 + 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 @@ -236,56 +240,51 @@ def read_skymap(params, map_struct=None): 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 + if 'do_3d' not in params: + 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) @@ -321,7 +320,7 @@ 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: ( @@ -333,20 +332,22 @@ 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_hdu: + 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 return params, map_struct From 3d53431fbd0e22460c92a4e23f49745d5961a8b4 Mon Sep 17 00:00:00 2001 From: Theodlz Date: Mon, 2 Dec 2024 15:29:57 -0800 Subject: [PATCH 2/3] single quotes to double quotes for consistency --- gwemopt/io/skymap.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gwemopt/io/skymap.py b/gwemopt/io/skymap.py index 8ae31b9..9f211f4 100644 --- a/gwemopt/io/skymap.py +++ b/gwemopt/io/skymap.py @@ -258,8 +258,8 @@ def read_skymap(params, map_struct=None, no_hdu=False): if params["gpstime"] is None: params["gpstime"] = t_obs.gps - if 'do_3d' not in params: - params['do_3d'] = do_3d + if "do_3d" not in params: + params["do_3d"] = do_3d map_struct = {} From 2d748452bfb6651f07cf25765664905cc7eeb61b Mon Sep 17 00:00:00 2001 From: Theodlz Date: Mon, 2 Dec 2024 16:01:36 -0800 Subject: [PATCH 3/3] remove no_hdu parameter, and replace by no_plotting, no_catalog, and no_efficiency parameters --- gwemopt/io/skymap.py | 46 +++++++++++++++++++++++++++++++++----------- 1 file changed, 35 insertions(+), 11 deletions(-) diff --git a/gwemopt/io/skymap.py b/gwemopt/io/skymap.py index 9f211f4..f326cc1 100644 --- a/gwemopt/io/skymap.py +++ b/gwemopt/io/skymap.py @@ -214,7 +214,7 @@ def read_inclination(skymap, params, map_struct): return map_struct -def read_skymap(params, map_struct=None, no_hdu=False): +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 @@ -224,6 +224,14 @@ def read_skymap(params, map_struct=None, no_hdu=False): :return: params, map_struct """ + # 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": @@ -292,15 +300,19 @@ def read_skymap(params, map_struct=None, no_hdu=False): 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() @@ -322,7 +334,14 @@ def read_skymap(params, map_struct=None, no_hdu=False): map_struct["skymap_raster_schedule"]["PROB"][ind[ii]] = 0.0 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"], @@ -333,7 +352,8 @@ def read_skymap(params, map_struct=None, no_hdu=False): ) map_struct["hdu"] = None - if not no_hdu: + 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"), @@ -350,4 +370,8 @@ def read_skymap(params, map_struct=None, no_hdu=False): 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