Skip to content

Commit

Permalink
mocpy 0.14.0
Browse files Browse the repository at this point in the history
  • Loading branch information
mcoughlin committed Jun 3, 2024
1 parent 6bfc1fd commit ef870e8
Show file tree
Hide file tree
Showing 8 changed files with 62 additions and 247 deletions.
23 changes: 13 additions & 10 deletions gwemopt/chipgaps/decam.py
Original file line number Diff line number Diff line change
Expand Up @@ -212,21 +212,24 @@ def get_decam_ccds(ra, dec, save_footprint=False):


def get_decam_quadrant_moc(ra, dec):
ccd_coords = get_decam_ccds(0, 0)
# ccd_coords = get_decam_ccds(0, 0)

skyoffset_frames = SkyCoord(ra, dec, unit=u.deg).skyoffset_frame()
ccd_coords_icrs = SkyCoord(
*np.tile(ccd_coords[:, np.newaxis, ...], (1, 1, 1)),
unit=u.deg,
frame=skyoffset_frames[:, np.newaxis, np.newaxis],
).transform_to(ICRS)
# skyoffset_frames = SkyCoord(ra, dec, unit=u.deg).skyoffset_frame()
# ccd_coords_icrs = SkyCoord(
# *np.tile(ccd_coords[:, np.newaxis, ...], (1, 1, 1)),
# unit=u.deg,
# frame=skyoffset_frames[:, np.newaxis, np.newaxis],
# ).transform_to(ICRS)

ccd_coords = get_decam_ccds(ra, dec)
ras, decs = ccd_coords[0], ccd_coords[1]

moc = None
for subfield_id, coords in enumerate(ccd_coords_icrs[0]):
for ra, dec in zip(ras, decs):
if moc is None:
moc = MOC.from_polygon_skycoord(coords)
moc = MOC.from_polygon(ra * u.deg, dec * u.deg)
else:
moc = moc + MOC.from_polygon_skycoord(coords)
moc = moc + MOC.from_polygon(ra * u.deg, dec * u.deg)

return moc

Expand Down
35 changes: 12 additions & 23 deletions gwemopt/coverage.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,11 @@

import ligo.segments as segments
import numpy as np
import regions
from astropy import units as u
from astropy.coordinates import SkyCoord, get_sun
from astropy.time import Time
from mocpy import MOC
from tqdm import tqdm

import gwemopt.plotting
Expand All @@ -22,7 +24,6 @@
slice_map_tiles,
slice_number_tiles,
)
from gwemopt.utils.pixels import getCirclePixels, getSquarePixels
from gwemopt.utils.treasuremap import get_treasuremap_pointings


Expand Down Expand Up @@ -87,32 +88,20 @@ def read_coverage(params, telescope, filename, moc_struct=None):
coverage_struct["filters"].append(filt)

if moc_struct is None:
if telescope == "ATLAS":
alpha = 0.2
color = "#6c71c4"
elif telescope == "PS1":
alpha = 0.1
color = "#859900"
else:
alpha = 0.2
color = "#6c71c4"

if config_struct["FOV_coverage_type"] == "square":
moc = getSquarePixels(
ra,
dec,
config_struct["FOV_coverage"],
alpha=alpha,
color=color,
center = SkyCoord(ra, dec, unit="deg", frame="icrs")
region = regions.RectangleSkyRegion(
center,
config_struct["FOV_coverage"] * u.deg,
config_struct["FOV_coverage"] * u.deg,
)
moc = MOC.from_astropy_regions(region, max_depth=10)
elif config_struct["FOV_coverage_type"] == "circle":
moc = getCirclePixels(
ra,
dec,
config_struct["FOV_coverage"],
alpha=alpha,
color=color,
center = SkyCoord(ra, dec, unit="deg", frame="icrs")
region = regions.CircleSkyRegion(
center, radius=config_struct["FOV"] * u.deg
)
moc = MOC.from_astropy_regions(region, max_depth=10)
else:
moc = moc_struct[field]["moc"]

Expand Down
27 changes: 14 additions & 13 deletions gwemopt/moc.py
Original file line number Diff line number Diff line change
@@ -1,17 +1,20 @@
import copy
import os

import astropy.units as u
import healpy as hp
import numpy as np
import regions
from astropy.coordinates import SkyCoord
from joblib import Parallel, delayed
from regions import Regions
from mocpy import MOC
from tqdm import tqdm

import gwemopt.tiles
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 getCirclePixels, getRegionPixels, getSquarePixels
from gwemopt.utils.pixels import getRegionPixels


def create_moc(params, map_struct=None):
Expand Down Expand Up @@ -114,24 +117,22 @@ def Fov2Moc(params, config_struct, telescope, ra_pointing, dec_pointing, nside):
moc_struct = {}

if config_struct["FOV_type"] == "square":
moc = getSquarePixels(
ra_pointing,
dec_pointing,
config_struct["FOV"],
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":
moc = getCirclePixels(
ra_pointing,
dec_pointing,
config_struct["FOV"],
)
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"])
regions = Regions.read(region_file, format="ds9")
region = regions.Regions.read(region_file, format="ds9")
moc = getRegionPixels(
ra_pointing,
dec_pointing,
regions,
region,
nside,
)
else:
Expand Down
26 changes: 8 additions & 18 deletions gwemopt/tiles.py
Original file line number Diff line number Diff line change
Expand Up @@ -1087,25 +1087,15 @@ def compute_tiles_map(
vals = np.nan * np.ones((ntiles,))

if params["doParallel"]:

def calculate_probability(moc, skymap):
return moc.probability_in_multiordermap(skymap)

with tqdm_joblib(
tqdm(desc="Probability calcuation", total=len(keys))
) as progress_bar:
vals = Parallel(
n_jobs=params["Ncores"],
backend=params["parallelBackend"],
batch_size=int(len(keys) / params["Ncores"]) + 1,
)(
delayed(calculate_probability)(tile_struct[key]["moc"], skymap)
for key in keys
)
vals = np.array(vals)
vals = MOC.probabilities_in_multiordermap(
[tile_struct[key]["moc"] for key in keys],
skymap,
n_threads=params["Ncores"],
)
else:
for ii, key in tqdm(enumerate(keys), total=len(keys)):
vals[ii] = tile_struct[key]["moc"].probability_in_multiordermap(skymap)
vals = MOC.probabilities_in_multiordermap(
[tile_struct[key]["moc"] for key in keys], skymap
)

return vals

Expand Down
7 changes: 1 addition & 6 deletions gwemopt/utils/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,6 @@
from gwemopt.utils.misc import auto_rasplit, get_exposures, integrationTime
from gwemopt.utils.observability import calculate_observability
from gwemopt.utils.param_utils import readParamsFromFile
from gwemopt.utils.pixels import (
get_ellipse_coords,
getCirclePixels,
getRectanglePixels,
getSquarePixels,
)
from gwemopt.utils.pixels import getRegionPixels
from gwemopt.utils.sidereal_time import greenwich_sidereal_time
from gwemopt.utils.treasuremap import get_treasuremap_pointings
172 changes: 0 additions & 172 deletions gwemopt/utils/pixels.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,34 +7,6 @@
from mocpy import MOC


def get_ellipse_coords(a=0.0, b=0.0, x=0.0, y=0.0, angle=0.0, npts=10):
"""Draws an ellipse using (360*k + 1) discrete points; based on pseudo code
given at http://en.wikipedia.org/wiki/Ellipse
k = 1 means 361 points (degree by degree)
a = major axis distance,
b = minor axis distance,
x = offset along the x-axis
y = offset along the y-axis
angle = clockwise rotation [in degrees] of the ellipse;
* angle=0 : the ellipse is aligned with the positive x-axis
* angle=30 : rotated 30 degrees clockwise from positive x-axis
"""
pts = np.zeros((npts, 2))

beta = -angle * np.pi / 180.0
sin_beta = np.sin(beta)
cos_beta = np.cos(beta)
alpha = np.linspace(0, 2 * np.pi, npts)

sin_alpha = np.sin(alpha)
cos_alpha = np.cos(alpha)

pts[:, 0] = x + (a * cos_alpha * cos_beta - b * sin_alpha * sin_beta)
pts[:, 1] = y + (a * cos_alpha * sin_beta + b * sin_alpha * cos_beta)

return pts


def getRegionPixels(
ra_pointing,
dec_pointing,
Expand Down Expand Up @@ -74,147 +46,3 @@ def getRegionPixels(
moc = moc + MOC.from_polygon_skycoord(coords)

return moc


def getCirclePixels(
ra_pointing,
dec_pointing,
radius,
):

radecs = get_ellipse_coords(
a=radius / np.cos(np.deg2rad(dec_pointing)),
b=radius,
x=ra_pointing,
y=dec_pointing,
angle=0.0,
npts=25,
)
idx = np.where(radecs[:, 1] > 90.0)[0]
radecs[idx, 1] = 180.0 - radecs[idx, 1]
idx = np.where(radecs[:, 1] < -90.0)[0]
radecs[idx, 1] = -180.0 - radecs[idx, 1]
idx = np.where(radecs[:, 0] > 360.0)[0]
radecs[idx, 0] = 720.0 - radecs[idx, 0]
idx = np.where(radecs[:, 0] < 0.0)[0]
radecs[idx, 0] = 360.0 + radecs[idx, 0]

radecs = np.array(radecs)
coords = coordinates.SkyCoord(radecs[:, 0] * u.deg, radecs[:, 1] * u.deg)
moc = MOC.from_polygon_skycoord(coords)

return moc


def getRectanglePixels(
ra_pointing,
dec_pointing,
raSide,
decSide,
):
area = raSide * decSide

decCorners = (dec_pointing - decSide / 2.0, dec_pointing + decSide / 2.0)

# security for the periodic limit conditions
radecs = []
for d in decCorners:
if d > 90.0:
d = 180.0 - d
elif d < -90.0:
d = -180 - d

raCorners = (
ra_pointing - (raSide / 2.0) / np.cos(np.deg2rad(d)),
ra_pointing + (raSide / 2.0) / np.cos(np.deg2rad(d)),
)
# security for the periodic limit conditions
for r in raCorners:
if r > 360.0:
r = r - 360.0
elif r < 0.0:
r = 360.0 + r
radecs.append([r, d])

radecs = np.array(radecs)
idx1 = np.where(radecs[:, 0] >= 180.0)[0]
idx2 = np.where(radecs[:, 0] < 180.0)[0]
idx3 = np.where(radecs[:, 0] > 300.0)[0]
idx4 = np.where(radecs[:, 0] < 60.0)[0]
if (len(idx1) > 0 and len(idx2) > 0) and not (len(idx3) > 0 and len(idx4) > 0):
alpha = 0.0

idx1 = np.where(np.abs(radecs[:, 1]) >= 87.0)[0]
if len(idx1) == 4:
return MOC.new_empty(max_depth=10)

idx1 = np.where((radecs[:, 1] >= 87.0) | (radecs[:, 1] <= -87.0))[0]
if len(idx1) > 0:
radecs = np.delete(radecs, idx1[0], 0)

npts, junk = radecs.shape
if npts == 4:
idx = [0, 1, 3, 2]
radecs = radecs[idx, :]

coords = coordinates.SkyCoord(radecs[:, 0] * u.deg, radecs[:, 1] * u.deg)
moc = MOC.from_polygon_skycoord(coords)

return moc


def getSquarePixels(
ra_pointing,
dec_pointing,
tileSide,
):
area = tileSide * tileSide

decCorners = (dec_pointing - tileSide / 2.0, dec_pointing + tileSide / 2.0)

# security for the periodic limit conditions
radecs = []
for d in decCorners:
if d > 90.0:
d = 180.0 - d
elif d < -90.0:
d = -180 - d

raCorners = (
ra_pointing - (tileSide / 2.0) / np.cos(np.deg2rad(d)),
ra_pointing + (tileSide / 2.0) / np.cos(np.deg2rad(d)),
)

# security for the periodic limit conditions
for r in raCorners:
if r > 360.0:
r = r - 360.0
elif r < 0.0:
r = 360.0 + r
radecs.append([r, d])

radecs = np.array(radecs)
idx1 = np.where(radecs[:, 0] >= 180.0)[0]
idx2 = np.where(radecs[:, 0] < 180.0)[0]
idx3 = np.where(radecs[:, 0] > 300.0)[0]
idx4 = np.where(radecs[:, 0] < 60.0)[0]
if (len(idx1) > 0 and len(idx2) > 0) and not (len(idx3) > 0 and len(idx4) > 0):
alpha = 0.0

idx1 = np.where(np.abs(radecs[:, 1]) >= 87.0)[0]
if len(idx1) == 4:
return MOC.new_empty(max_depth=10)

idx1 = np.where((radecs[:, 1] >= 87.0) | (radecs[:, 1] <= -87.0))[0]
if len(idx1) > 0:
radecs = np.delete(radecs, idx1[0], 0)

npts, junk = radecs.shape
if npts == 4:
idx = [0, 1, 3, 2]
radecs = radecs[idx, :]

coords = coordinates.SkyCoord(radecs[:, 0] * u.deg, radecs[:, 1] * u.deg)
moc = MOC.from_polygon_skycoord(coords)

return moc
Loading

0 comments on commit ef870e8

Please sign in to comment.