Merge pull request #772 from lincc-frameworks/from_catalog

Improve documentation for position sampling

Author: jeremykubica

docs/faq.rst

@@ -16,7 +16,7 @@ your simulation:
 You can also use `logging.DEBUG` to see even more detailed debug-level messages.
 
 
-I Am Seeing Empty Light Curves for Some Entrys. Why Is This Happening?
+I Am Seeing Empty Light Curves for Some Entries. Why Is This Happening?
 -------------------------------------------------------------------------------
 
 The most common reason that LightCurveLynx produces rows with empty light curves is
@@ -80,3 +80,16 @@ Yes with some caveats. LightCurveLynx has built-in support for simulating spectr
 in an early stage of development and does not yet add noise to the measurements. In addition, spectra
 simulation is **only** compaible with models that generate data on the spectral level (not bandflux-only
 models).  For more detail see :doc:`the spectrograph demo notebook <notebooks/spectrograph_demo>`.
+
+
+Can I Generate Points from a Catalog?
+--------------------------------------------------------------------------------
+
+Yes. LightCurveLynx allows you to generate light curves for objects in a catalog containing 
+positions using the``CatalogRADECSampler`` object This sampler takes in a table
+of information with at least "ra" and "dec" columns. The helper function 
+``from_hats()`` is provided to load directly from a [HATS](https://www.ivoa.net/documents/Notes/HATS/)
+catalog.
+
+See the :doc:`sampling positions demo notebook <notebooks/sampling_positions>` notebook
+for a detailed description of how to sample (RA, dec) positions.

docs/notebooks/sampling_positions.ipynb

@@ -6,11 +6,16 @@
    "source": [
     "# Sampling (RA, dec)\n",
     "\n",
-    "LightCurveLynx provides multiple mechanisms for sampling (RA, dec) from choosing points uniformly on a sphere to using the footprint of a given survey. In this notebook we discuss several of the approaches and their relative tradeoffs.\n",
+    "Sampling the position (RA, dec) of an object on the sky is a critical for in survey-based simulation. LightCurveLynx matches the sampled positions (RA, dec) with the survey information (`ObsTable`) to determine when each object is observed. Points from the object's light curve are *only* generated at those times where the object is observed. For example if an object is at (45.0, -10.0) and the survey includes this points at MJDs 60676.0 and 60678.0, the generated light curve will only have two fluxes. More importantly, if the sampled position falls outside the survey's footprint the returned light curve will be empty (the object is never observed). Therefore it is critical to choose a reasonable sampling scheme.\n",
     "\n",
-    "LightCurveLynx matches the sampled positions (RA, dec) with the survey information (`ObsTable`) to determine when each object is observed. Points from the object's light curve are *only* generated at those times where the object is observed. For example if an object is at (45.0, -10.0) and the survey includes this points at MJDs 60676.0 and 60678.0, the generated light curve will only have two fluxes. More importantly, if the sampled position falls outside the survey's footprint the returned light curve will be empty (the object is never observed). Therefor it is critical to choose a reasonable sampling scheme.\n",
+    "LightCurveLynx provides multiple mechanisms for sampling object positions, including:\n",
     "\n",
-    "For almost all cases we recommend the `ApproximateMOCSampler` for generating positions that correspond to the underlying survey."
+    "  * Complete uniform over a sphere (`UniformRADEC`)\n",
+    "  * Sampling from a survey (`ObsTableRADECSampler`, `ObsTableUniformRADECSampler`\n",
+    ", and `ApproximateMOCSampler`)\n",
+    "  * Sampling from a catalog (`CatalogRADECSampler`)\n",
+    "\n",
+    "The choice of sampler will depend heavily on the science use case."
    ]
   },
   {
@@ -21,9 +26,11 @@
    "source": [
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
+    "import pandas as pd\n",
     "\n",
     "from lightcurvelynx.math_nodes.ra_dec_sampler import (\n",
     "    ApproximateMOCSampler,\n",
+    "    CatalogRADECSampler,\n",
     "    ObsTableRADECSampler,\n",
     "    ObsTableUniformRADECSampler,\n",
     "    UniformRADEC,\n",
@@ -60,15 +67,17 @@
    "source": [
     "However this approach has limited use when simulating a specific survey. Depending on the survey's coverage, a significant number of (RA, dec) points may fall outside the viewing area. In those cases, the returned light curves will be empty, indicating that there were no observations of a given object.\n",
     "\n",
-    "## Sampling from a Survey\n",
+    "## Sampling from a Survey's Footprint\n",
     "\n",
     "We can sample (RA, dec) coordinates from a survey (an `ObsTable` object) in two ways. First we could sample a pointing from the survey and then a point from that field of view. Second, we could sample uniformly from the region coverage by the survey.\n",
     "\n",
-    "We consider each of these approaches below, but recommend the `ApproximateMOCSampler` for the majority of simulations that want to generate positions from a given footprint.\n",
+    "We consider each of these approaches below, but recommend the `ApproximateMOCSampler` for the majority of simulations that want to generate positions from a given footprint as it is most efficient over a range of scenarios.\n",
+    "\n",
+    "### Sampling Pointings (Visted Weighted)\n",
     "\n",
-    "### Sampling Pointings\n",
+    "Sampling pointings from the survey provides a **visit-weighted** sampling of positions covered by the survey. The sampling works from a table of all the survey's pointings and randomly selecting a row for each sample.\n",
     "\n",
-    "Sampling pointings from the survey provides a **visit-weighted** sampling of positions covered by the survey. For concreteness let's start with a survey that visits two fields: one centered at (45.0, -15.0) and the other at (315.0, 15.0). The first field is visited once and the second field is visited four times on four consecutive nights."
+    "For concreteness let's start with a survey that visits two fields: one centered at (45.0, -15.0) and the other at (315.0, 15.0). The first field is visited once and the second field is visited four times on four consecutive nights."
    ]
   },
   {
@@ -110,7 +119,10 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "As we can see, the field centered in the Northern hemisphere is sampled significantly more than the one centered in the Southern hemisphere.\n",
+    "As we can see, the field centered in the Northern hemisphere is sampled significantly more than the one centered in the Southern hemisphere, because the survey visited that field more often.\n",
+    "\n",
+    "Users can limit the overlap with the `dedup_threshold` parameter, which removes near duplicate rows. But for non-visit-weighted sampling of a survey, we highly recommend using one of either `ObsTableUniformRADECSampler` or `ApproximateMOCSampler`.\n",
+    "\n",
     "\n",
     "### Sampling Survey Coverage\n",
     "\n",
@@ -169,6 +181,9 @@
    "source": [
     "As you can see, many of the points land outside the 1 degree radius around the center of the pointing (10, 0).\n",
     "\n",
+    "We recommend the `ObsTableUniformRADECSampler` **only** for cases where there is large survey coverage.\n",
+    "\n",
+    "\n",
     "**ApproximateMOCSampler**\n",
     "\n",
     "The `ApproximateMOCSampler` samples from the area covered by a [Multi-Order Coverage Map (MOC)](https://www.ivoa.net/documents/MOC/20190215/WD-MOC-1.1-20190215.pdf), which is a collection of healpix pixels representing an area on the sky. Users can generate custom MOCs for hypothetical surveys, build a MOC from a survey, or use a helper function to create the sampler directly from the survey. **This is our recommended approach for sampling from a survey.**\n",
@@ -219,7 +234,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "For example, if we use depth=8, the survey's coverage is approximated by a grid of only 786,432 pixels over the entire sky with (an average pixel width of around 14 arc minutes). We recommend at least a depth of 12 (average pixel width around 50 arc seconds) for reasonable accuracy.\n",
+    "For example, if we use `depth=8`, the survey's coverage is approximated by a grid of only 786,432 pixels over the entire sky with (an average pixel width of around 14 arc minutes). We recommend at least a depth of 12 (average pixel width around 50 arc seconds) for reasonable accuracy.\n",
     "\n",
     "As a concrete example, let's look at what happens if we prebuild the MOC at depth=4 (very coarse). Although the sampler uses a depth of 14, it cannot extract any more resolution from the input than it was given (a depth=4 MOC)."
    ]
@@ -275,6 +290,44 @@
     "moc_sampler.plot_footprint()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Sampling from a Catalog\n",
+    "\n",
+    "We can use the `CatalogRADECSampler` node to sample from a catalog of known objects. This class is a thin connivence wrapper on the `ObsTableRADECSampler` that samples only ra and dec and defaults to a zero radius (exact sampling of observations). It provides a `dedup_threshold` parameter for removing near duplicates."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 4 points with one pair of duplicates\n",
+    "values = {\n",
+    "    \"ra\": np.array([45.0, 315.0, 315.0, 0.0]),\n",
+    "    \"dec\": np.array([-15.0, 15.0, 15.0, -10.0]),\n",
+    "}\n",
+    "data = pd.DataFrame(values)\n",
+    "\n",
+    "pointing_sampler = CatalogRADECSampler(data, dedup_threshold=0.1, node_label=\"catalog\")\n",
+    "print(f\"The sampler will draw from {len(pointing_sampler)} points.\")\n",
+    "\n",
+    "(ra, dec) = pointing_sampler.generate(num_samples=100)\n",
+    "\n",
+    "plt.scatter(ra, dec, s=1)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can use the `CatalogRADECSampler` (as well as the `ObsTableRADECSampler`) includes a `from_hats()` function that allows users to directly load a catalog from the [HATS format](https://www.ivoa.net/documents/Notes/HATS/)."
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

src/lightcurvelynx/math_nodes/ra_dec_sampler.py

@@ -4,11 +4,13 @@
 from pathlib import Path
 
 import numpy as np
+import pandas as pd
 from astropy.coordinates import Angle, SkyCoord
 from cdshealpix.nested import healpix_to_skycoord
-from citation_compass import CiteClass
+from citation_compass import CiteClass, cite_inline
 from mocpy import MOC
 
+from lightcurvelynx.astro_utils.coordinate_utils import dedup_coords
 from lightcurvelynx.math_nodes.given_sampler import TableSampler
 from lightcurvelynx.math_nodes.np_random import NumpyRandomFunc
 from lightcurvelynx.obstable.obs_table import ObsTable
@@ -96,11 +98,15 @@ class ObsTableRADECSampler(TableSampler):
         Use 0.0 to return the exact given points.
         If None and data is an ObsTable, uses the value from the ObsTable.
         Default: None
+    dedup_threshold : float, optional
+        The deduplication threshold in degrees. If two rows have RA and dec values that
+        are within this threshold, only one of them will be kept for sampling. Use 0.0 to
+        keep all rows. Default: 0.0
     **kwargs : dict, optional
         Additional keyword arguments to pass to the parent class constructor.
     """
 
-    def __init__(self, data, *, extra_cols=None, radius=None, **kwargs):
+    def __init__(self, data, *, extra_cols=None, radius=None, dedup_threshold=0.0, **kwargs):
         if isinstance(data, ObsTable):
             if radius is None:
                 radius = data.radius
@@ -112,6 +118,9 @@ def __init__(self, data, *, extra_cols=None, radius=None, **kwargs):
             raise ValueError(f"Invalid radius: {radius}")
         self.radius = radius
 
+        if "ra" not in data or "dec" not in data:
+            raise ValueError("Data must contain 'ra' and 'dec' columns.")
+
         # Start with RA, dec, and (optionally) time.
         data_dict = {
             "ra": data["ra"],
@@ -126,11 +135,29 @@ def __init__(self, data, *, extra_cols=None, radius=None, **kwargs):
                 if col not in data_dict:
                     data_dict[col] = data[col]
 
+        # Do a deduplication step to remove rows with very close RA and dec values.
+        if dedup_threshold > 0.0:
+            data_dict = pd.DataFrame(data_dict)
+            _, _, inds = dedup_coords(
+                data_dict["ra"].values,
+                data_dict["dec"].values,
+                threshold=dedup_threshold,
+            )
+            data_dict = data_dict.iloc[inds].reset_index(drop=True)
+
         super().__init__(data_dict, in_order=False, **kwargs)
 
     @classmethod
-    def from_hats(cls, path, *, radius=None, extra_cols=None, **kwargs):
-        """Create a GivenRADECSampler from the observations in a HATS Catalog.
+    def from_hats(
+        cls,
+        path,
+        *,
+        radius=None,
+        extra_cols=None,
+        dedup_threshold=0.0,
+        **kwargs,
+    ):
+        """Create a ObsTableRADECSampler from the observations in a HATS Catalog.
 
         Note
         ----
@@ -148,13 +175,17 @@ def from_hats(cls, path, *, radius=None, extra_cols=None, **kwargs):
         extra_cols : list of str, optional
             A list of extra column names to include in the sampling.
             Default: None
+        dedup_threshold : float, optional
+            The deduplication threshold in degrees. If two rows have RA and dec values that
+            are within this threshold, only one of them will be kept for sampling. Use 0.0 to
+            keep all rows. Default: 0.0
         **kwargs : dict, optional
             Additional keyword arguments to pass to the constructor.
 
         Returns
         -------
-        GivenRADECSampler
-            The created GivenRADECSampler object.
+        ObsTableRADECSampler
+            The created ObsTableRADECSampler object.
         """
         # See if the (optional) LSDB package is installed.
         try:
@@ -171,8 +202,17 @@ def from_hats(cls, path, *, radius=None, extra_cols=None, **kwargs):
             cols_to_load.extend(extra_cols)
         columns = list(set(cols_to_load))  # Remove any duplicates.
 
+        # Cite the HATS format.
+        cite_inline("HATS Catalog Format", "https://www.ivoa.net/documents/Notes/HATS/)")
+
         data = read_hats(path, columns=columns).compute()
-        return cls(data, extra_cols=extra_cols, radius=radius, **kwargs)
+        return cls(
+            data,
+            extra_cols=extra_cols,
+            radius=radius,
+            dedup_threshold=dedup_threshold,
+            **kwargs,
+        )
 
     def compute(self, graph_state, rng_info=None, **kwargs):
         """Return the given values.
@@ -660,3 +700,29 @@ def compute(self, graph_state, rng_info=None, **kwargs):
         graph_state.set(self.node_string, "ra", ra)
         graph_state.set(self.node_string, "dec", dec)
         return (ra, dec)
+
+
+class CatalogRADECSampler(ObsTableRADECSampler):
+    """A FunctionNode that randomly samples RA and dec from a given catalog of objects.
+
+    Note
+    ----
+    This is a thin convenience wrapper around ObsTableRADECSampler.
+
+    Parameters
+    ----------
+    data : Pandas DataFrame, NestedFrame, or dict
+        The data to use for sampling. Must contain 'ra' and 'dec' columns.
+    dedup_threshold : float, optional
+        The deduplication threshold in degrees. If two rows have RA and dec values that
+        are within this threshold, only one of them will be kept for sampling. Use 0.0 to
+        keep all rows. Default: 0.0
+    **kwargs : dict, optional
+        Additional keyword arguments to pass to the parent class constructor.
+    """
+
+    def __init__(self, data, *, dedup_threshold=0.0, **kwargs):
+        # Always default to a radius of 0.0.
+        if "radius" not in kwargs or kwargs["radius"] is None:
+            kwargs["radius"] = 0.0
+        super().__init__(data, dedup_threshold=dedup_threshold, **kwargs)