Skip to content

refactor: dedupe FootprintPlotter, import from cs_util #186

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
cailmdaley wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

refactor: dedupe FootprintPlotter, import from cs_util #186

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion pyproject.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ dependencies = [
"camb",
"clmm",
"colorama",
"cs_util>=0.1.5",
"cs_util>=0.1.9",
"emcee",
"getdist @ git+https://github.com/benabed/getdist.git@upper_triangle_whisker",
"h5py",
Expand Down
273 changes: 1 addition & 272 deletions src/sp_validation/plots.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
from astropy import units as u
from astropy.coordinates import SkyCoord
from cs_util import plots
from cs_util.plots import FootprintPlotter
from lenspack.geometry.projections.gnom import radec2xy

from sp_validation.plot_style import *
Expand Down Expand Up @@ -399,278 +400,6 @@ def plot_binned(
plots.savefig(f"{key}_binned.png")


class FootprintPlotter:
"""Class to create footprint plots.

Parameters
-----------
nside_coverage: int, optional
basic resolution of map; default is 32
nside_map:
fine resolution for plotting; default is 2048

"""

# Dictionary storing region parameters
_regions = {
"NGC": {"ra_0": 180, "extend": [120, 270, 20, 70], "vmax": 60},
"SGC": {"ra_0": 15, "extend": [-20, 45, 20, 45], "vmax": 60},
"fullsky": {"ra_0": 150, "extend": [0, 360, -90, 90], "vmax": 60},
}

def __init__(self, nside_coverage=32, nside_map=2048):

self._nside_coverage = nside_coverage
self._nside_map = nside_map

def create_hsp_map(self, ra, dec):
"""Create Hsp Map.

Create healsparse map.

Parameters
----------
ra : numpy.ndarray
right ascension values
dec : numpy.ndarray
declination values

Returns
-------
hsp.HealSparseMap
map

"""
# Create empty map
hsp_map = hsp.HealSparseMap.make_empty(
self._nside_coverage, self._nside_map, dtype=np.float32, sentinel=np.nan
)

# Get pixel list corresponding to coordinates
hpix = hp.ang2pix(self._nside_map, ra, dec, nest=True, lonlat=True)

# Get count of objects per pixel
pixel_counts = Counter(hpix)

# List of unique pixels
unique_hpix = np.array(list(pixel_counts.keys()))

# Number of objects
values = np.array(list(pixel_counts.values()), dtype=np.float32)

# Create maps with numbers per pixel
hsp_map[unique_hpix] = values

return hsp_map

def plot_area(
self,
hsp_map,
ra_0=0,
extend=[120, 270, 29, 70],
vmax=60,
projection=None,
outpath=None,
title=None,
colorbar=True,
colorbar_label="Coverage depth",
):
"""Plot Area.

Plot catalogue in an area on the sky.

Parameters
----------
hsp_map : hsp_HealSparseMap
input map
ra_0 : float, optional
anchor point in R.A.; default is 0
extend : list, optional
sky region, extend=[ra_low, ra_high, dec_low, dec_high];
default is [120, 270, 29, 70]
vmax : float, optional
maximum pixel value to plot with color; default is 60
projection : skyproj.McBrydeSkyproj
if ``None`` (default), a new plot is created
outpath : str, optional
output path, default is ``None``
title : str, optional
print title if not ``None`` (default)
colorbar : bool, optional
add colorbar; default is ``True``
colorbar_label : str, optional
colorbar label; default is "Coverage depth"

Returns
--------
skyproj.McBrydeSkyproj
projection instance
plt.axes.Axes
axes instance

Raises
------
ValueError
if no object found in region

"""
if not projection:

# Create new figure and axes
fig, ax = plt.subplots(figsize=(10, 10))

# Create new projection
projection = skyproj.McBrydeSkyproj(
ax=ax, lon_0=ra_0, extent=extend, autorescale=True, vmax=vmax
)
else:
ax = None

im = None
try:
im, lon_raster, lat_raster, values_raster = projection.draw_hspmap(
hsp_map, lon_range=extend[0:2], lat_range=extend[2:]
)
except ValueError:
msg = "No object found in region to draw"
print(f"{msg}, continuing...")
# raise ValueError(msg)

projection.draw_milky_way(width=25, linewidth=1.5, color="black", linestyle="-")

# Add colorbar if requested and image was drawn
if colorbar and im is not None:
plt.colorbar(im, ax=ax if ax else projection.ax, label=colorbar_label)

if title:
plt.title(title, pad=5)

if outpath:
plt.savefig(outpath)

return projection, ax

def plot_region(self, hsp_map, region, projection=None, outpath=None, title=None, colorbar=True, colorbar_label="Coverage depth"):
"""Plot Region.

Plot catalogue in a predefined region on the sky.

Parameters
----------
hsp_map : hsp_HealSparseMap
input map
region : dict
region dictionary with keys 'ra_0', 'extend', 'vmax'
projection : skyproj.McBrydeSkyproj, optional
if ``None`` (default), a new plot is created
outpath : str, optional
output path, default is ``None``
title : str, optional
print title if not ``None`` (default)
colorbar : bool, optional
add colorbar; default is ``True``
colorbar_label : str, optional
colorbar label; default is "Coverage depth"

Returns
--------
skyproj.McBrydeSkyproj
projection instance
plt.axes.Axes
axes instance

"""
return self.plot_area(
hsp_map,
region["ra_0"],
region["extend"],
region["vmax"],
projection=projection,
outpath=outpath,
title=title,
colorbar=colorbar,
colorbar_label=colorbar_label,
)

def plot_all_regions(self, hsp_map, outbase=None):

for region in self._regions:
if outbase:
outpath = f"{outbase}_{region}.png"
else:
outpath = None
self.plot_region(hsp_map, self._regions[region], outpath=outpath)

@classmethod
def hp_pixel_centers(cls, nside, nest=False):

# Get number of pixels for given nside
npix = hp.nside2npix(nside)

# Get pixel indices
pix_indices = np.arange(npix)

# Get coordinates of pixel centers
ra, dec = hp.pix2ang(nside, pix_indices, nest=nest, lonlat=True)

return ra, dec, npix

@classmethod
def plot_footprint_as_hp(cls, hsp_map, nside, outpath=None, title=None):

ra, dec, npix = cls.hp_pixel_centers(nside)

# Create an empty HEALPix map
m = np.full(npix, np.nan)

fig, ax = plt.subplots(figsize=(10, 10))

# Plot the HEALPix grid
hp.mollview(m, title=title, coord="C", notext=True, rot=(180, 0, 0))

# Define the Galactic Plane: l = [0, 360], b = 0°
for l0, ls in zip((-5, 0, 5), (":", "-", ":")):
l_values = np.linspace(0, 360, 500) # 500 points along the plane
b_values = np.zeros_like(l_values) # Galactic latitude is 0 (the plane)

# Convert (l, b) to (λ, β) - Ecliptic coordinates
coords = SkyCoord(
l=l_values * u.degree, b=b_values * u.degree, frame="galactic"
)
ecl_coords = coords.transform_to(
"barycentrictrueecliptic"
) # Ecliptic frame

# Extract Ecliptic longitude (λ) and latitude (β)
lambda_ecl = ecl_coords.lon.deg # Ecliptic longitude
beta_ecl = ecl_coords.lat.deg # Ecliptic latitude

# Convert to HEALPix projection coordinates (colatitude, longitude)
theta = np.radians(90 - beta_ecl) # HEALPix uses colatitude
phi = np.radians(lambda_ecl) # HEALPix uses longitude

# Create a healpy Mollweide projection in Ecliptic coordinates
hp.projplot(
theta, phi, linestyle=ls, color="black", linewidth=1
) # Plot the outline

# Apply mask
mask_values = hsp_map.get_values_pos(ra, dec, valid_mask=True, lonlat=True)

ok = np.where(mask_values == False)[0]
# nok = np.where(mask_values == False)[0]

hp.projscatter(ra[ok], dec[ok], lonlat=True, color="green", s=1, marker=".")
# hp.projscatter(ra[nok], dec[nok], lonlat=True, color="red", s=1, marker=".")

plt.tight_layout()

if outpath:
plt.savefig(outpath)

plt.show()


def hsp_map_logical_or(maps, verbose=False):
"""
Hsp Map Logical Or.
Expand Down
Loading