From 288baa540add3063b83c2f9dfce69bef0be90885 Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Tue, 25 Jun 2024 04:49:01 -0500 Subject: [PATCH] MOC.from_cones and MOC.from_boxes --- gwemopt/moc.py | 138 +++++++++++--------------------------- gwemopt/utils/__init__.py | 2 +- gwemopt/utils/pixels.py | 58 ++++++---------- 3 files changed, 62 insertions(+), 136 deletions(-) diff --git a/gwemopt/moc.py b/gwemopt/moc.py index c135538..852331b 100644 --- a/gwemopt/moc.py +++ b/gwemopt/moc.py @@ -14,7 +14,7 @@ from gwemopt.chipgaps import get_decam_quadrant_moc, get_ztf_quadrant_moc from gwemopt.paths import CONFIG_DIR from gwemopt.utils.parallel import tqdm_joblib -from gwemopt.utils.pixels import getRegionPixels +from gwemopt.utils.pixels import get_region_moc def create_moc(params, map_struct=None): @@ -26,6 +26,13 @@ def create_moc(params, map_struct=None): tesselation = config_struct["tesselation"] moc_struct = {} + if (telescope == "ZTF") and params["doUsePrimary"]: + idx = np.where(tesselation[:, 0] <= 880)[0] + tesselation = tesselation[idx, :] + elif (telescope == "ZTF") and params["doUseSecondary"]: + idx = np.where(tesselation[:, 0] >= 1000)[0] + tesselation = tesselation[idx, :] + if params["doChipGaps"]: if params["doParallel"]: n_threads = params["Ncores"] @@ -42,63 +49,41 @@ def create_moc(params, map_struct=None): else: raise ValueError("Chip gaps only available for DECam and ZTF") - for ii, tess in enumerate(tesselation): - index, ra, dec = tess[0], tess[1], tess[2] - if (telescope == "ZTF") and params["doUsePrimary"] and (index > 880): - continue - if (telescope == "ZTF") and params["doUseSecondary"] and (index < 1000): - continue - index = index.astype(int) - moc_struct[index] = {"ra": ra, "dec": dec, "moc": mocs[ii]} else: if params["doParallel"]: - with tqdm_joblib( - tqdm(desc="MOC creation", total=len(tesselation)) - ) as progress_bar: - moclists = Parallel( - n_jobs=params["Ncores"], - backend=params["parallelBackend"], - batch_size=int(len(tesselation) / params["Ncores"]) + 1, - )( - delayed(Fov2Moc)( - params, config_struct, telescope, tess[1], tess[2], nside - ) - for tess in tesselation - ) - for ii, tess in tqdm(enumerate(tesselation), total=len(tesselation)): - index, ra, dec = tess[0], tess[1], tess[2] - if ( - (telescope == "ZTF") - and params["doUsePrimary"] - and (index > 880) - ): - continue - if ( - (telescope == "ZTF") - and params["doUseSecondary"] - and (index < 1000) - ): - continue - moc_struct[index] = moclists[ii] + n_threads = params["Ncores"] else: - for ii, tess in tqdm(enumerate(tesselation), total=len(tesselation)): - index, ra, dec = tess[0], tess[1], tess[2] - if ( - (telescope == "ZTF") - and params["doUsePrimary"] - and (index > 880) - ): - continue - if ( - (telescope == "ZTF") - and params["doUseSecondary"] - and (index < 1000) - ): - continue - index = index.astype(int) - moc_struct[index] = Fov2Moc( - params, config_struct, telescope, ra, dec, nside - ) + n_threads = None + + if config_struct["FOV_type"] == "circle": + mocs = MOC.from_cones( + lon=tesselation[:, 1] * u.deg, + lat=tesselation[:, 2] * u.deg, + radius=config_struct["FOV"] * u.deg, + max_depth=np.uint8(10), + n_threads=n_threads, + ) + elif config_struct["FOV_type"] == "square": + mocs = MOC.from_boxes( + lon=tesselation[:, 1] * u.deg, + lat=tesselation[:, 2] * u.deg, + a=config_struct["FOV"] * u.deg, + b=config_struct["FOV"] * u.deg, + angle=0 * u.deg, + max_depth=np.uint8(10), + n_threads=n_threads, + ) + elif config_struct["FOV_type"] == "region": + region_file = os.path.join(CONFIG_DIR, config_struct["FOV"]) + region = regions.Regions.read(region_file, format="ds9") + mocs = get_region_moc( + tesselation[:, 1], tesselation[:, 2], region, n_threads=n_threads + ) + + moc_struct = { + tess[0].astype(int): {"ra": tess[1], "dec": tess[2], "moc": mocs[ii]} + for ii, tess in enumerate(tesselation) + } if map_struct is not None: moc_keep = map_struct["moc_keep"] @@ -141,46 +126,3 @@ def create_moc(params, map_struct=None): moc_structs[telescope] = moc_struct return moc_structs - - -def Fov2Moc(params, config_struct, telescope, ra_pointing, dec_pointing, nside): - """Return a MOC in fits file of a fov footprint. - The MOC fov is displayed in real time in an Aladin plan. - - Input: - ra--> right ascention of fov center [deg] - dec --> declination of fov center [deg] - fov_width --> fov width [deg] - fov_height --> fov height [deg] - nside --> healpix resolution; by default - """ - - moc_struct = {} - - if config_struct["FOV_type"] == "square": - center = SkyCoord(ra_pointing, dec_pointing, unit="deg", frame="icrs") - region = regions.RectangleSkyRegion( - center, config_struct["FOV"] * u.deg, config_struct["FOV"] * u.deg - ) - moc = MOC.from_astropy_regions(region, max_depth=10) - elif config_struct["FOV_type"] == "circle": - center = SkyCoord(ra_pointing, dec_pointing, unit="deg", frame="icrs") - region = regions.CircleSkyRegion(center, radius=config_struct["FOV"] * u.deg) - moc = MOC.from_astropy_regions(region, max_depth=10) - elif config_struct["FOV_type"] == "region": - region_file = os.path.join(CONFIG_DIR, config_struct["FOV"]) - region = regions.Regions.read(region_file, format="ds9") - moc = getRegionPixels( - ra_pointing, - dec_pointing, - region, - nside, - ) - else: - raise ValueError("FOV_type must be square, circle or region") - - moc_struct["ra"] = ra_pointing - moc_struct["dec"] = dec_pointing - moc_struct["moc"] = moc - - return moc_struct diff --git a/gwemopt/utils/__init__.py b/gwemopt/utils/__init__.py index 1998752..e7aba99 100644 --- a/gwemopt/utils/__init__.py +++ b/gwemopt/utils/__init__.py @@ -3,6 +3,6 @@ from gwemopt.utils.misc import get_exposures, integrationTime from gwemopt.utils.observability import calculate_observability from gwemopt.utils.param_utils import readParamsFromFile -from gwemopt.utils.pixels import getRegionPixels +from gwemopt.utils.pixels import get_region_moc from gwemopt.utils.sidereal_time import greenwich_sidereal_time from gwemopt.utils.treasuremap import get_treasuremap_pointings diff --git a/gwemopt/utils/pixels.py b/gwemopt/utils/pixels.py index 411ca57..40dedbb 100644 --- a/gwemopt/utils/pixels.py +++ b/gwemopt/utils/pixels.py @@ -1,48 +1,32 @@ import astropy.units as u -import healpy as hp -import matplotlib import numpy as np -from astropy import coordinates -from astropy_healpix import HEALPix -from mocpy import MOC +from astropy.coordinates import ICRS, SkyCoord +from mocpy import MOC, mocpy +from tqdm import tqdm -def getRegionPixels( - ra_pointing, - dec_pointing, - regions, - nside, -): - theta = 0.5 * np.pi - np.deg2rad(dec_pointing) - phi = np.deg2rad(ra_pointing) +def get_region_moc(ra, dec, regions, max_depth=12, n_threads=None): - HPX = HEALPix(nside=nside, order="nested", frame=coordinates.ICRS()) - - skyoffset_frames = coordinates.SkyCoord( - ra_pointing, dec_pointing, unit=u.deg - ).skyoffset_frame() - - moc = None - - # security for the periodic limit conditions for reg in regions: ra_tmp = reg.vertices.ra dec_tmp = reg.vertices.dec coords = np.stack([np.array(ra_tmp), np.array(dec_tmp)]) - # Copy the tile coordinates such that there is one per field - # in the grid - coords_icrs = coordinates.SkyCoord( - *np.tile(coords[:, np.newaxis, ...], (1, 1, 1)), - unit=u.deg, - frame=skyoffset_frames[:, np.newaxis, np.newaxis], - ).transform_to(coordinates.ICRS) - coords = coords_icrs.transform_to(HPX.frame) - - if moc is None: - moc = MOC.from_polygon_skycoord(coords) - else: - moc = moc + MOC.from_polygon_skycoord(coords) - - return moc + skyoffset_frames = SkyCoord(ra, dec, unit=u.deg).skyoffset_frame() + coords_icrs = SkyCoord( + *np.tile(coords[:, np.newaxis, ...], (1, 1, 1)), + unit=u.deg, + frame=skyoffset_frames[:, np.newaxis, np.newaxis], + ).transform_to(ICRS) + + mocs = [] + for ccd_coords in tqdm(coords_icrs): + stacked = np.stack((ccd_coords.ra.deg, ccd_coords.dec.deg), axis=1) + result = stacked.reshape(-1, ccd_coords.ra.deg.shape[1]) + lon_lat_list = [row for row in result] + indices = mocpy.from_polygons(lon_lat_list, np.uint8(max_depth), n_threads) + moc = sum([MOC(index) for index in indices]) + mocs.append(moc) + + return mocs