Skip to content

Commit

Permalink
MOC.from_cones and MOC.from_boxes
Browse files Browse the repository at this point in the history
  • Loading branch information
@mcoughlin
mcoughlin committed Jun 25, 2024
1 parent b3f2e39 commit 288baa5
Show file tree
Hide file tree
Showing 3 changed files with 62 additions and 136 deletions.
138 changes: 40 additions & 98 deletions gwemopt/moc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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"]
Expand All @@ -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"]
Expand Down Expand Up @@ -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
2 changes: 1 addition & 1 deletion gwemopt/utils/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
58 changes: 21 additions & 37 deletions gwemopt/utils/pixels.py
View file
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit 288baa5

Please sign in to comment.