diff --git a/gwemopt/io/skymap.py b/gwemopt/io/skymap.py index bb39444..f326cc1 100644 --- a/gwemopt/io/skymap.py +++ b/gwemopt/io/skymap.py @@ -214,21 +214,33 @@ 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 @@ -236,20 +248,18 @@ 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 @@ -257,35 +267,32 @@ def read_skymap(params, map_struct=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) @@ -293,15 +300,19 @@ def read_skymap(params, map_struct=None): 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() @@ -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"], @@ -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