feat: Add coordinate_system parameter to compute_tangential_and_cross_components and _compute_lensing_angles functions

clmm/dataops/__init__.py

@@ -1,4 +1,5 @@
 """Functions to compute polar/azimuthal averages in radial bins"""
+
 import warnings
 import numpy as np
 import scipy
@@ -14,6 +15,7 @@
     _validate_ra,
     _validate_dec,
     _validate_is_deltasigma_sigma_c,
+    _validate_coordinate_system,
 )
 from ..redshift import (
     _integ_pzfuncs,
@@ -29,6 +31,7 @@ def compute_tangential_and_cross_components(
     dec_source,
     shear1,
     shear2,
+    coordinate_system="pixel",
     geometry="curve",
     is_deltasigma=False,
     sigma_c=None,
@@ -95,6 +98,9 @@ def compute_tangential_and_cross_components(
         The measured shear (or reduced shear or ellipticity) of the source galaxies
     shear2: array
         The measured shear (or reduced shear or ellipticity) of the source galaxies
+    coordinate_system: str, optional
+        Coordinate system of the ellipticity components. Options are 'pixel' or 'sky'.
+        Default is 'pixel'.
     geometry: str, optional
         Sky geometry to compute angular separation.
         Options are curve (uses astropy) or flat.
@@ -128,6 +134,7 @@ def compute_tangential_and_cross_components(
         validate_argument(locals(), "shear2", "float_array")
         validate_argument(locals(), "geometry", str)
         validate_argument(locals(), "sigma_c", "float_array", none_ok=True)
+        _validate_coordinate_system(locals(), "coordinate_system", str)
         ra_source_, dec_source_, shear1_, shear2_ = arguments_consistency(
             [ra_source, dec_source, shear1, shear2],
             names=("Ra", "Dec", "Shear1", "Shear2"),
@@ -147,9 +154,13 @@ def compute_tangential_and_cross_components(
         )
     # Compute the lensing angles
     if geometry == "flat":
-        angsep, phi = _compute_lensing_angles_flatsky(ra_lens, dec_lens, ra_source_, dec_source_)
+        angsep, phi = _compute_lensing_angles_flatsky(
+            ra_lens, dec_lens, ra_source_, dec_source_, coordinate_system=coordinate_system
+        )
     elif geometry == "curve":
-        angsep, phi = _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_, dec_source_)
+        angsep, phi = _compute_lensing_angles_astropy(
+            ra_lens, dec_lens, ra_source_, dec_source_, coordinate_system=coordinate_system
+        )
     else:
         raise NotImplementedError(f"Sky geometry {geometry} is not currently supported")
     # Compute the tangential and cross shears
@@ -329,7 +340,9 @@ def compute_galaxy_weights(
     return w_ls
 
 
-def _compute_lensing_angles_flatsky(ra_lens, dec_lens, ra_source_list, dec_source_list):
+def _compute_lensing_angles_flatsky(
+    ra_lens, dec_lens, ra_source_list, dec_source_list, coordinate_system="pixel"
+):
     r"""Compute the angular separation between the lens and the source and the azimuthal
     angle from the lens to the source in radians.
 
@@ -353,6 +366,9 @@ def _compute_lensing_angles_flatsky(ra_lens, dec_lens, ra_source_list, dec_sourc
         Right ascensions of each source galaxy in degrees
     dec_source_list: array
         Declinations of each source galaxy in degrees
+    coordinate_system: str, optional
+        Coordinate system of the ellipticity components. Options are 'pixel' or 'sky'.
+        Default is 'pixel'.
 
     Returns
     -------
@@ -376,12 +392,16 @@ def _compute_lensing_angles_flatsky(ra_lens, dec_lens, ra_source_list, dec_sourc
         phi[angsep == 0.0] = 0.0
     else:
         phi = 0.0 if angsep == 0.0 else phi
+    if coordinate_system == "sky":
+        phi = np.pi - phi
     if np.any(angsep > np.pi / 180.0):
         warnings.warn("Using the flat-sky approximation with separations >1 deg may be inaccurate")
     return angsep, phi
 
 
-def _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_list, dec_source_list):
+def _compute_lensing_angles_astropy(
+    ra_lens, dec_lens, ra_source_list, dec_source_list, coordinate_system="pixel"
+):
     r"""Compute the angular separation between the lens and the source and the azimuthal
     angle from the lens to the source in radians.
 
@@ -414,6 +434,8 @@ def _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_list, dec_sourc
     else:
         phi -= 2 * np.pi if phi > np.pi else 0
         phi = 0 if angsep == 0 else phi
+    if coordinate_system == "sky":
+        phi = np.pi - phi
     return angsep, phi
 
 

clmm/galaxycluster.py

@@ -1,6 +1,7 @@
 """@file galaxycluster.py
 The GalaxyCluster class
 """
+
 import pickle
 import warnings
 from .gcdata import GCData
@@ -17,6 +18,7 @@
     _validate_ra,
     _validate_dec,
     _draw_random_points_from_tab_distribution,
+    _validate_coordinate_system,
 )
 
 
@@ -35,6 +37,9 @@ class GalaxyCluster:
         Redshift of galaxy cluster center
     galcat : GCData
         Table of background galaxy data containing at least galaxy_id, ra, dec, e1, e2, z
+    coordinate_system : str
+        Coordinate system of the galaxy cluster center (pixel or sky)
+
     validate_input: bool
         Validade each input argument
     """
@@ -51,13 +56,22 @@ def __init__(self, *args, validate_input=True, **kwargs):
             self._check_types()
             self.set_ra_lower(ra_low=0)
 
-    def _add_values(self, unique_id: str, ra: float, dec: float, z: float, galcat: GCData):
+    def _add_values(
+        self,
+        unique_id: str,
+        ra: float,
+        dec: float,
+        z: float,
+        galcat: GCData,
+        coordinate_system: str = "pixel",
+    ):
         """Add values for all attributes"""
         self.unique_id = unique_id
         self.ra = ra
         self.dec = dec
         self.z = z
         self.galcat = galcat
+        self.coordinate_system = coordinate_system
 
     def _check_types(self):
         """Check types of all attributes"""
@@ -66,6 +80,7 @@ def _check_types(self):
         _validate_dec(vars(self), "dec", False)
         validate_argument(vars(self), "z", (float, str), argmin=0, eqmin=True)
         validate_argument(vars(self), "galcat", GCData)
+        _validate_coordinate_system(vars(self), "coordinate_system", str)
         self.unique_id = str(self.unique_id)
         self.ra = float(self.ra)
         self.dec = float(self.dec)
@@ -242,6 +257,9 @@ def compute_tangential_and_cross_components(
             Name of the column to be added to the `galcat` astropy table that will contain the
             cross component computed from columns `shape_component1` and `shape_component2`.
             Default: `ex`
+        coordinate_system: str, optional
+            Coordinate system of the ellipticity components. Options are 'pixel' or 'sky'.
+            Default: 'pixel'
         geometry: str, optional
             Sky geometry to compute angular separation.
             Options are curve (uses astropy) or flat.
@@ -281,6 +299,7 @@ def compute_tangential_and_cross_components(
             dec_lens=self.dec,
             geometry=geometry,
             validate_input=self.validate_input,
+            coordinate_system=self.coordinate_system,
             **cols,
         )
         if add:

clmm/support/mock_data.py

@@ -38,6 +38,7 @@ def generate_galaxy_catalog(
     ngal_density=None,
     pz_bins=101,
     pzpdf_type="shared_bins",
+    coordinate_system="pixel",
     validate_input=True,
 ):
     r"""Generates a mock dataset of sheared background galaxies.
@@ -156,6 +157,9 @@ def generate_galaxy_catalog(
         The number density of galaxies (in galaxies per square arcminute, from z=0 to z=infty).
         The number of galaxies to be drawn will then depend on the redshift distribution and
         user-defined redshift range.  If specified, the ngals argument will be ignored.
+    coordinate_system : str, optional
+        The coordinate system to use for the output catalog. Options are 'pixel' and 'sky'.
+
     validate_input: bool
         Validade each input argument
 
@@ -231,6 +235,7 @@ def generate_galaxy_catalog(
         "pz_bins": pz_bins,
         "field_size": field_size,
         "pzpdf_type": pzpdf_type,
+        "coordinate_system": coordinate_system,
     }
 
     if ngals is None and ngal_density is None:
@@ -361,6 +366,7 @@ def _generate_galaxy_catalog(
     photoz_sigma_unscaled=None,
     pz_bins=101,
     pzpdf_type="shared_bins",
+    coordinate_system="pixel",
     field_size=None,
 ):
     """A private function that skips the sanity checks on derived properties. This
@@ -422,6 +428,9 @@ def _generate_galaxy_catalog(
     _, posangle = c_cl.separation(c_gal).rad, c_cl.position_angle(c_gal).rad
     posangle += 0.5 * np.pi  # for right convention
 
+    if coordinate_system == "sky":
+        posangle = np.pi - posangle  # ellipticity coordinate system conversion
+
     # corresponding shear1,2 components
     gam1 = -gamt * np.cos(2 * posangle) + gamx * np.sin(2 * posangle)
     gam2 = -gamt * np.sin(2 * posangle) - gamx * np.cos(2 * posangle)

clmm/utils/__init__.py

@@ -40,6 +40,7 @@
     _validate_ra,
     _validate_dec,
     _validate_is_deltasigma_sigma_c,
+    _validate_coordinate_system,
 )
 
 from .units import (

clmm/utils/validation.py

@@ -1,4 +1,5 @@
 """General utility functions that are used in multiple modules"""
+
 import numpy as np
 from ..constants import Constants as const
 
@@ -224,3 +225,23 @@ def _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c):
         raise TypeError("sigma_c (=None) must be provided when is_deltasigma=True")
     if not is_deltasigma and sigma_c is not None:
         raise TypeError(f"sigma_c (={sigma_c}) must be None when is_deltasigma=False")
+
+
+def _validate_coordinate_system(loc, coordinate_system, valid_type):
+    r"""Validate the coordinate system.
+
+    Parameters
+    ----------
+    loc: dict
+        Dictionary with all input arguments. Should be locals().
+    coordinate_system: str
+        Coordinate system used for the input data.
+    valid_type: str, type
+        Valid types for argument, options are object types, list/tuple of types, or:
+
+            * 'int_array' - interger, interger array
+            * 'float_array' - float, float array
+    """
+    validate_argument(loc, coordinate_system, valid_type)
+    if loc[coordinate_system] not in ["sky", "pixel"]:
+        raise ValueError(f"{coordinate_system} must be 'sky' or 'pixel'.")