Skip to content

Commit

Permalink
Ruff and docstring updates
Browse files Browse the repository at this point in the history
  • Loading branch information
rhiannonlynne committed Jan 15, 2024
1 parent 2d07f04 commit 68d271e
Show file tree
Hide file tree
Showing 17 changed files with 331 additions and 247 deletions.
8 changes: 4 additions & 4 deletions rubin_sim/maf/maps/base_map.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ def help(cls, doc=False):


class BaseMap(metaclass=MapsRegistry):
""" """
"""Base for maps classes."""

def __init__(self, **kwargs):
self.keynames = ["newkey"]
Expand All @@ -65,8 +65,8 @@ def __ge__(self, othermap):
return self.keynames >= othermap.keynames

def run(self, slice_points):
"""
Given slice_points (dict containing metadata about each slice_point, including ra/dec),
adds additional metadata at each slice_point and returns updated dict.
"""Given slice_points (dict containing metadata about each slice_point,
including ra/dec), adds additional metadata at each slice_point
and returns updated dict.
"""
raise NotImplementedError("This must be defined in subclass")
33 changes: 21 additions & 12 deletions rubin_sim/maf/maps/create_gaia_density_map.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,17 +7,21 @@

# Modifying createStarDensityMap to use GAIA DR1 catalog

# Use the catsim framework to loop over a healpy map and generate a stellar density map
# Use the catsim framework to loop over a healpy map and generate a
# stellar density map

# Connect to fatboy with: ssh -L 51433:fatboy.phys.washington.edu:1433 gateway.astro.washington.edu
# Connect to fatboy with: ssh -L 51433:fatboy.phys.washington.edu:1433
# gateway.astro.washington.edu
# If non-astro user, use [email protected]

# NOTE: fatboy is no longer operative

if __name__ == "__main__":
# Hide imports here so documentation builds
from rubin_sim.catalogs.db import DBObject
from rubin_sim.utils import angular_separation, halfSpaceFromRaDec
from lsst.sims.catalogs.db import DBObject
from lsst.sims.utils import angular_separation, halfSpaceFromRaDec

# from rubin_sim.catalogs.generation.db import CatalogDBObject
# from lsst.sims.catalogs.generation.db import CatalogDBObject
# Import the bits needed to get the catalog to work
# from rubin_sim.catUtils.baseCatalogModels import *
# from rubin_sim.catUtils.exampleCatalogDefinitions import *
Expand Down Expand Up @@ -61,7 +65,7 @@
over_max_mask = data["over_max_mask"].copy()

print("")
# Look at a cirular area the same area as the healpix it's centered on.
# Look at a circular area the same area as the healpix it's centered on.
bound_length = hpsize_deg / np.pi**0.5
radius = bound_length

Expand All @@ -76,12 +80,14 @@
chunk_size = 10000
for i in np.arange(indx_min, int(npix)):
last_cp = ""
# wonder what the units of bound_length are...degrees! And it's a radius
# wonder what the units of bound_length are...degrees!
# And it's a radius
# The newer interface:
# obs_metadata = ObservationMetaData(bound_type='circle',
## pointing_ra=np.degrees(ra[i]),
# pointing_dec=np.degrees(dec[i]),
# bound_length=bound_length, mjd=5700)
# bound_length=bound_length,
# mjd=5700)

# t = dbobj.getCatalog('ref_catalog_star', obs_metadata=obs_metadata)
hs = halfSpaceFromRaDec(ra[i], dec[i], radius)
Expand All @@ -106,18 +112,21 @@
results = gaia_db.get_arbitrary_chunk_iterator(query, dtype=dtype, chunk_size=10000)
result = list(results)[0]

distances = angular_separation(result["ra"], result["dec"], ra[i], dec[i]) # Degrees
distances = angular_separation(result["ra"], result["dec"], ra[i], dec[i])
result = result[np.where(distances < radius)]

import pdb

pdb.set_trace()
# I could think of setting the chunksize to something really large, then only doing one chunk?
# Or maybe setting up a way to break out of the loop if everything gets really dense?
# I could think of setting the chunksize to something really large,
# then only doing one chunk?
# Or maybe setting up a way to break out of the loop if
# everything gets really dense?
temp_hist = np.zeros(np.size(bins) - 1, dtype=float)
counter = 0
col_name = "phot_g_mean_mag"
for chunk in results:
chunk_hist, bins = np.histogram(chunk[colName], bins)
chunk_hist, bins = np.histogram(chunk[col_name], bins)
temp_hist += chunk_hist
counter += chunk_size
if counter >= break_limit:
Expand Down
27 changes: 15 additions & 12 deletions rubin_sim/maf/maps/dust_map.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,37 +8,39 @@


class DustMap(BaseMap):
"""
Compute the E(B-V) for each point in a given spatial distribution of slicePoints.
"""Add the E(B-V) values to the slice points.
Primarily, this calls eb_vhp to read a healpix map of E(B-V) values over
the sky, then assigns ebv values to each slice_point.
If the slicer is a healpix slicer, this is trivial.
Otherwise, it either uses the nearest healpix grid point or interpolates.
Primarily, this calls eb_vhp to read a healpix map of E(B-V) values over the sky, then
assigns ebv values to each slice_point. If the slicer is a healpix slicer, this is trivial.
The key added to the slice points is `ebv`.
Parameters
----------
interp : `bool`, opt
Interpolate the dust map at each slice_point (True) or just use the nearest value (False).
Interpolate the dust map at each slice_point (True)
or just use the nearest value (False).
Default is False.
nside : `int`, opt
Default nside value to read the dust map from disk. Primarily useful if the slicer is not
a healpix slicer.
Default nside value to read the dust map from disk.
Primarily useful if the slicer is not a healpix slicer.
Default 128.
map_path : `str`, opt
Define a path to the directory holding the dust map files.
Default None, which uses RUBIN_SIM_DATA_DIR.
"""

def __init__(self, interp=False, nside=128, map_path=None):
"""
interp: should the dust map be interpolated (True) or just use the nearest value (False).
"""
self.keynames = ["ebv"]
self.interp = interp
self.nside = nside
self.map_path = map_path

def run(self, slice_points):
# If the slicer has nside, it's a healpix slicer so we can read the map directly
# If the slicer has nside,
# it's a healpix slicer so we can read the map directly
if "nside" in slice_points:
if slice_points["nside"] != self.nside:
warnings.warn(
Expand All @@ -50,7 +52,8 @@ def run(self, slice_points):
pixels=slice_points["sid"],
map_path=self.map_path,
)
# Not a healpix slicer, look up values based on RA,dec with possible interpolation
# Not a healpix slicer,
# look up values based on RA,dec with possible interpolation
else:
slice_points["ebv"] = eb_vhp(
self.nside,
Expand Down
84 changes: 52 additions & 32 deletions rubin_sim/maf/maps/dust_map_3d.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,38 +11,48 @@


class DustMap3D(BaseMap):
"""The DustMap3D provides a `~rubin_sim.maf.map` to hold 3d EBV data.
Adds the following keys to the slicePoints:
ebv3d_dists - the distances from the 3d dust map at each slice_point (in pc)
ebv3d_ebvs - the E(B-V) values corresponding to each distance at each slice_point
ebv3d_ebv_at_<dist_pc> - the (single) ebv value at the nearest distance to dist_pc
ebv3d_dist_at_<d_mag> - the (single) distance value corresponding to where extinction and
distance modulus combine to create a m-Mo value of d_mag, for the filter specified in filtername (in pc).
Note that <dist_pc> and <d_mag> will be formatted with a single decimal place.
The additional method 'distance_at_mag' can be called either with the distances and ebv values for the
entire map or with the values from a single slice_point, in order to calculate the distance at which
extinction and distance modulus combine to create a m-Mo value closest to 'dmag' in any filter.
This is the same value as would be reported in ebv3d_dist_at_<d_mag>, but can be calculated on the fly,
"""Add 3-d E(B-V) values to the slice points.
The slice point dictionary keys are expanded with the following keys:
ebv3d_dists -
the distances from the 3d dust map at each slice_point (in pc)
`ebv3d_ebvs` -
the E(B-V) values corresponding to each distance at each slice_point
`ebv3d_ebv_at_<dist_pc>` -
the (single) ebv value at the nearest distance to dist_pc
`ebv3d_dist_at_<d_mag>` -
the (single) distance value corresponding to where extinction and
distance modulus combine to create a m-Mo value of d_mag, for the filter
specified in filtername (in pc).
Note that <dist_pc> and <d_mag> will be formatted with a
single decimal place.
The additional method 'distance_at_mag' can be called either with the
distances and ebv values for the entire map or with the values from a
single slice_point, in order to calculate the distance at which
extinction and distance modulus combine to create a m-Mo value closest
to 'dmag' in any filter. This is the same value as would be reported in
ebv3d_dist_at_<d_mag>, but can be calculated on the fly,
allowing variable filters and dmag values.
Parameters
----------
nside: `int`
Healpixel resolution (2^x).
nside : `int`
Healpixel resolution (2^x) to read from disk.
map_file : `str`, opt
Path to dust map file.
interp : `bool`, opt
Should returned values be interpolated (True) or just nearest neighbor (False).
Should returned values be interpolated (True)
or just nearest neighbor (False).
Default True, but is ignored if 'pixels' is provided.
filtername : 'str', opt
Name of the filter (to match the lsst filter names in rubin_sim.photUtils.DustValues)
in which to calculate dust extinction magnitudes
Name of the filter (to match the lsst filter names in
rubin_sim.photUtils.DustValues) in which to calculate dust
extinction magnitudes
dist_pc : `float`, opt
Distance at which to precalculate the nearest ebv value
Distance at which to precalculate the nearest ebv value (pc)
d_mag : `float`, opt
Calculate the maximum distance which matches this 'd_mag'
Calculate the maximum distance which matches this `d_mag`
d_mag == m-mO (dust extinction + distance modulus)
r_x : `dict` {`str`: `float`}, opt
Per-filter dust extinction curve coefficients.
Expand All @@ -65,10 +75,12 @@ def __init__(
self.filtername = filtername
self.dist_pc = dist_pc
self.d_mag = d_mag
# r_x is the extinction coefficient (A_v = R_v * E(B-V) .. A_x = r_x * E(B-V)) per filter
# This is equivalent to calculating A_x (using rubin_sim.photUtils.Sed.addDust) in each
# filter and setting E(B-V) to 1 [so similar to the values calculated in DustValues ..
# we probably should rename those (from Ax1 to r_x)
# r_x is the extinction coefficient (A_v = R_v * E(B-V) ..
# A_x = r_x * E(B-V)) per filter
# This is equivalent to calculating A_x
# (using rubin_sim.photUtils.Sed.addDust) in each
# filter and setting E(B-V) to 1 [so similar to the values
# calculated in DustValues.
if r_x is None:
self.r_x = DustValues().r_x.copy()
else:
Expand All @@ -82,7 +94,8 @@ def __init__(
]

def run(self, slice_points):
# If the slicer has nside, it's a healpix slicer so we can read the map directly
# If the slicer has nside,
# it's a healpix slicer so we can read the map directly
if "nside" in slice_points:
if slice_points["nside"] != self.nside:
warnings.warn(
Expand All @@ -94,7 +107,8 @@ def run(self, slice_points):
pixels=slice_points["sid"],
map_file=self.map_file,
)
# Not a healpix slicer, look up values based on RA,dec with possible interpolation
# Not a healpix slicer,
# look up values based on RA,dec with possible interpolation
else:
dists, ebvs = ebv_3d_hp(
self.nside,
Expand All @@ -107,7 +121,8 @@ def run(self, slice_points):
# Calculate the map ebv and dist values at the initialized distance
dist_closest, ebv_at_dist = get_x_at_nearest_y(dists, ebvs, self.dist_pc)

# Calculate the distances at which m_minus_Mo values of 'dmag' are reached
# Calculate the distances at which m_minus_Mo values
# of 'dmag' are reached
dist_dmag = self.distance_at_dmag(self.d_mag, dists, ebvs, self.filtername)

slice_points["ebv3d_dists"] = dists
Expand All @@ -118,21 +133,26 @@ def run(self, slice_points):
return slice_points

def distance_at_dmag(self, dmag, dists, ebvs, filtername=None):
# Provide this as a method which could be used for a single slice_point as well as for whole map
# (single slice_point means you could calculate this for any arbitrary magnitude or filter if needed)
# Provide this as a method which could be used for a single
# slice_point as well as for whole map
# (single slice_point means you could calculate this for any
# arbitrary magnitude or filter if needed)
# This method is here because some metrics require it.
if filtername is None:
filtername = self.filtername
# calculate distance modulus for each distance
dmods = 5.0 * np.log10(dists) - 5.0
# calculate dust extinction at each distance, for the filtername
a_x = self.r_x[filtername] * ebvs
# calculate the (m-Mo) = distmod + a_x -- combination of extinction due to distance and dust
# calculate the (m-Mo) = distmod + a_x -- combination of extinction
# due to distance and dust
m_minus__mo = dmods + a_x

# Maximum distance for the given m-Mo (dmag) value
# first do the 'within the map' closest distance
m_minus__mo_at_mag, dist_closest = get_x_at_nearest_y(m_minus__mo, dists, dmag)
# calculate distance modulus for an object with the maximum dust extinction (and then the distance)
# calculate distance modulus for an object with the maximum dust
# extinction (and then the distance)
if a_x.ndim == 2:
dist_mods_far = dmag - a_x[:, -1]
else:
Expand Down
31 changes: 20 additions & 11 deletions rubin_sim/maf/maps/ebv_3d_hp.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,25 +19,31 @@ def ebv_3d_hp(
pixels=None,
interp=False,
):
"""Reads and saves a 3d dust extinction file, return extinction at specified points (ra/dec/ or pixels).
"""Reads and saves a 3d dust extinction file from disk, return extinction
at specified points (ra/dec/ or pixels).
Parameters
----------
nside: `int`
nside : `int`
Healpixel resolution (2^x).
map_file : `str`, opt
Path to dust map file.
ra : `np.ndarray` or `float`, opt
RA (can take numpy array). Default None sets up healpix array of nside. Radians.
RA (can take numpy array).
Default None sets up healpix array of nside. Radians.
dec : `np.ndarray` or `float`, opt
Dec (can take numpy array). Default None set up healpix array of nside. Radians.
Dec (can take numpy array).
Default None set up healpix array of nside. Radians.
pixels : `np.ndarray`, opt
Healpixel IDs, to sub-select particular healpix points. Default uses all points.
Healpixel IDs, to sub-select particular healpix points.
Default uses all points.
Easiest way to access healpix values.
Note that the pixels in the healpix array MUST come from a healpix grid with the same nside
as the ebv_3d_hp map. Using different nsides can potentially fail silently.
Note that the pixels in the healpix array MUST come from a h
ealpix grid with the same nside as the ebv_3d_hp map.
Using different nsides can potentially fail silently.
interp : `bool`, opt
Should returned values be interpolated (True) or just nearest neighbor (False).
Should returned values be interpolated (True)
or just nearest neighbor (False).
Default False.
"""
if (ra is None) & (dec is None) & (pixels is None):
Expand Down Expand Up @@ -89,14 +95,16 @@ def ebv_3d_hp(
f"Will use nside from map data."
)
if pixels is not None:
# We're just going to raise an exception here because this could mean bad things.
# We're just going to raise an exception here
# because this could mean bad things.
raise ValueError(
f"Map nside {map_nside} did not match expected nside {nside}, "
f"and pixels provided; this can potentially indicate a serious "
f"error. Make nsides match or specify ra/dec instead of pixels."
)
nside = map_nside
# Nested healpix data will not match the healpix arrays for the slicers (at this time)
# Nested healpix data will not match the healpix arrays
# for the slicers (at this time)
if nested:
warnings.warn("Map has nested (not ring order) data; will reorder.")
for i in np.arange(0, len(dists[0])):
Expand Down Expand Up @@ -140,7 +148,8 @@ def get_x_at_nearest_y(x, y, x_goal):
the x at a single point of the map (1d array)
y : `np.array`
Can be either a map with y at each point in the map (2d array) or
the y at a single point of the map (1d array) - but should match x dimensionality
the y at a single point of the map (1d array) -
but should match x dimensionality
x_goal : `float'
The goal x value to look for the nearest y value
Expand Down
Loading

0 comments on commit 68d271e

Please sign in to comment.