Issue/635/Modify coordinate system implementation

clmm/dataops/__init__.py

@@ -1,4 +1,5 @@
 """Data operation for polar/azimuthal averages in radial bins and weights"""
+
 import warnings
 import numpy as np
 from astropy.coordinates import SkyCoord
@@ -25,7 +26,7 @@ def compute_tangential_and_cross_components(
     dec_source,
     shear1,
     shear2,
-    coordinate_system="euclidean",
+    coordinate_system=None,
     geometry="curve",
     is_deltasigma=False,
     sigma_c=None,
@@ -95,7 +96,7 @@ def compute_tangential_and_cross_components(
     coordinate_system: str, optional
         Coordinate system of the ellipticity components. Must be either 'celestial' or
         euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details.
-        Default is 'euclidean'.
+        If not set, defaults to 'euclidean'.
     geometry: str, optional
         Sky geometry to compute angular separation.
         Options are curve (uses astropy) or flat.
@@ -120,6 +121,10 @@ def compute_tangential_and_cross_components(
     # pylint: disable-msg=too-many-locals
     # Note: we make these quantities to be np.array so that a name is not passed from astropy
     # columns
+    if coordinate_system is None:
+        warnings.warn("The coordinate_system argument was not provided. Defaulting to 'euclidean'.")
+        coordinate_system = "euclidean"
+
     if validate_input:
         _validate_ra(locals(), "ra_source", True)
         _validate_dec(locals(), "dec_source", True)
@@ -129,7 +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)
+        _validate_coordinate_system(locals(), "coordinate_system")
         ra_source_, dec_source_, shear1_, shear2_ = arguments_consistency(
             [ra_source, dec_source, shear1, shear2],
             names=("Ra", "Dec", "Shear1", "Shear2"),
@@ -482,6 +487,7 @@ def make_radial_profile(
     z_lens=None,
     validate_input=True,
     weights=None,
+    coordinate_system=None,
 ):
     r"""Compute the angular profile of given components
 
@@ -535,6 +541,10 @@ def make_radial_profile(
     weights: array-like, optional
         Array of individual galaxy weights. If specified, the radial binned profile is
         computed using a weighted average
+    coordinate_system: str, optional
+        Coordinate system of the ellipticity components. Must be either 'celestial' or
+        euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details.
+        Default is 'euclidean'.
 
     Returns
     -------
@@ -553,6 +563,9 @@ def make_radial_profile(
     module.
     """
     # pylint: disable-msg=too-many-locals
+    if not coordinate_system:
+        warnings.warn("coordinate_system not set, defaulting to 'euclidean'.")
+        coordinate_system = "euclidean"
     if validate_input:
         validate_argument(locals(), "angsep", "float_array")
         validate_argument(locals(), "angsep_units", str)
@@ -564,6 +577,7 @@ def make_radial_profile(
         arguments_consistency(components, names=comp_dict.keys(), prefix="Input components")
         for component in comp_dict:
             validate_argument(comp_dict, component, "float_array")
+        _validate_coordinate_system(locals(), "coordinate_system")
     # Check to see if we need to do a unit conversion
     if angsep_units is not bin_units:
         source_seps = convert_units(angsep, angsep_units, bin_units, redshift=z_lens, cosmo=cosmo)
@@ -576,7 +590,7 @@ def make_radial_profile(
     profile_table = GCData(
         [bins[:-1], np.zeros(len(bins) - 1), bins[1:]],
         names=("radius_min", "radius", "radius_max"),
-        meta={"bin_units": bin_units},  # Add metadata
+        meta={"bin_units": bin_units, "coordinate_system": coordinate_system},  # Add metadata
     )
     # Compute the binned averages and associated errors
     for i, component in enumerate(components):

clmm/galaxycluster.py

@@ -18,7 +18,6 @@
     _validate_ra,
     _validate_dec,
     _draw_random_points_from_tab_distribution,
-    _validate_coordinate_system,
 )
 
 
@@ -37,10 +36,6 @@ 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, optional
-        Coordinate system of the ellipticity components. Must be either 'celestial' or
-        euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details.
-        Default is 'euclidean'.
     validate_input: bool
         Validade each input argument
     """
@@ -64,15 +59,13 @@ def _add_values(
         dec: float,
         z: float,
         galcat: GCData,
-        coordinate_system: str = "euclidean",
     ):
         """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"""
@@ -81,7 +74,6 @@ 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)
@@ -297,7 +289,7 @@ def compute_tangential_and_cross_components(
             dec_lens=self.dec,
             geometry=geometry,
             validate_input=self.validate_input,
-            coordinate_system=self.coordinate_system,
+            coordinate_system=self.galcat.meta["coordinate_system"],
             **cols,
         )
         if add:
@@ -607,6 +599,7 @@ def make_radial_profile(
                 for n in (tan_component_in_err, cross_component_in_err, None)
             ],
             weights=self.galcat[weights_in].data if use_weights else None,
+            coordinate_system=self.galcat.meta["coordinate_system"],
         )
         # Reaname table columns
         for i, name in enumerate([tan_component_out, cross_component_out, "z"]):
@@ -722,3 +715,20 @@ def set_ra_lower(self, ra_low=0):
         if "ra" in self.galcat.columns:
             self.galcat["ra"][self.galcat["ra"] < ra_low] += 360.0
             self.galcat["ra"][self.galcat["ra"] >= ra_low + 360.0] -= 360.0
+
+    def update_coordinate_system(self, coordinate_system, *args):
+        """Updates coordinate_system metadata of the galcat and converts ellipticity
+
+        Parameters
+        ----------
+        coordinate_system : str
+            Coordinate system of the ellipticity components. Must be either 'celestial' or
+            euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details.
+        *args : tuple
+            Components to be converted, must be in data.
+
+        Returns
+        -------
+        None
+        """
+        self.galcat.update_coordinate_system(coordinate_system, *args)

clmm/gcdata.py

@@ -1,16 +1,18 @@
 """
 Define the custom data type
 """
+
 import warnings
 from collections import OrderedDict
 from astropy.table import Table as APtable
-import numpy as np
-
 import qp
+import numpy as np
+from .utils import _validate_coordinate_system
 
 
 class GCMetaData(OrderedDict):
-    r"""Object to store metadata, it always has a cosmo key with protective changes
+    r"""Object to store metadata, it always has a cosmo and coordinate_system keys with protective
+    changes
 
     Attributes
     ----------
@@ -29,6 +31,10 @@ def __setitem__(self, item, value):
             raise ValueError(
                 "cosmo must be changed via update_cosmo or update_cosmo_ext_valid method"
             )
+        if item == "coordinate_system" and self.get("coordinate_system"):
+            raise ValueError(
+                "coordinate_system must be changed via update_coordinate_system method"
+            )
         OrderedDict.__setitem__(self, item, value)
 
     def __getitem__(self, item):
@@ -43,7 +49,7 @@ def __getitem__(self, item):
 
 class GCData(APtable):
     """
-    GCData: A data objetc for gcdata. Right now it behaves as an astropy table, with the following
+    GCData: A data object for gcdata. Right now it behaves as an astropy table, with the following
     modifications: `__getitem__` is case independent;
     The attribute .meta['cosmo'] is protected and
     can only be changed via update_cosmo or update_cosmo_ext_valid methods;
@@ -65,6 +71,16 @@ def __init__(self, *args, **kwargs):
         APtable.__init__(self, *args, **kwargs)
         metakwargs = kwargs["meta"] if "meta" in kwargs else {}
         metakwargs = {} if metakwargs is None else metakwargs
+        if (
+            "coordinate_system" not in metakwargs
+            and "copy_indices" not in kwargs
+            and "masked" not in kwargs
+            and "copy" not in kwargs
+        ):
+            warnings.warn("coordinate_system not set, defaulting to 'euclidean'")
+            metakwargs["coordinate_system"] = "euclidean"
+        if "coordinate_system" in metakwargs:
+            _validate_coordinate_system(metakwargs, "coordinate_system")
         self.meta = GCMetaData(**metakwargs)
         # this attribute is set when source galaxies have p(z)
         self.pzpdf_info = {
@@ -89,7 +105,7 @@ def _str_pzpdf_info(self):
                 out += " " + str(np.round(self.pzpdf_info.get("zbins"), 2))
                 np.set_printoptions(**default_cfg)
             elif out == "quantiles":
-                np.set_printoptions(formatter={'float': "{0:g}".format}, edgeitems=3, threshold=6)
+                np.set_printoptions(formatter={"float": "{0:g}".format}, edgeitems=3, threshold=6)
                 out += " " + str(self.pzpdf_info["quantiles"])
                 out += " - unpacked with zgrid : " + str(
                     self.pzpdf_info["unpack_quantile_zbins_limits"]
@@ -171,6 +187,8 @@ def update_info_ext_valid(self, key, gcdata, ext_value, overwrite=False):
         -------
         None
         """
+        if key == "coordinate_system":
+            _validate_coordinate_system(locals(), "ext_value")
         if ext_value:
             in_value = gcdata.meta[key]
             if in_value and in_value != ext_value:
@@ -221,6 +239,30 @@ def update_cosmo(self, cosmo, overwrite=False):
         """
         self.update_cosmo_ext_valid(self, cosmo, overwrite=overwrite)
 
+    def update_coordinate_system(self, coordinate_system, *args):
+        r"""Updates coordinate_system metadata and converts ellipticity
+
+        Parameters
+        ----------
+        coordinate_system : str
+            Coordinate system of the ellipticity components. Must be either 'celestial' or
+            euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details.
+        *args : tuple
+            Components to be converted, must be in data.
+
+        Returns
+        -------
+        None
+        """
+        in_coordinate_system = self.meta["coordinate_system"]
+        self.update_info_ext_valid("coordinate_system", self, coordinate_system, overwrite=True)
+        if coordinate_system != in_coordinate_system:
+            for col in args:
+                if col in self.colnames:
+                    self[col] = -self[col]
+                else:
+                    raise ValueError(f"component {col} not found in data")
+
     def has_pzpdfs(self):
         """Get pzbins and pzpdfs of galaxies
 

clmm/support/mock_data.py

@@ -44,7 +44,7 @@ def generate_galaxy_catalog(
     pz_bins=101,
     pz_quantiles_conf=(5, 31),
     pzpdf_type="shared_bins",
-    coordinate_system="euclidean",
+    coordinate_system=None,
     validate_input=True,
 ):
     r"""Generates a mock dataset of sheared background galaxies.
@@ -170,7 +170,7 @@ def generate_galaxy_catalog(
     coordinate_system : str, optional
         Coordinate system of the ellipticity components. Must be either 'celestial' or
         euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details.
-        Default is 'euclidean'.
+        If not set, defaults to 'euclidean'.
 
     validate_input: bool
         Validade each input argument
@@ -188,6 +188,10 @@ def generate_galaxy_catalog(
     # Too many local variables (25/15)
     # pylint: disable=R0914
 
+    if not coordinate_system:
+        warnings.warn("The coordinate_system argument was not provided. Defaulting to 'euclidean'.")
+        coordinate_system = "euclidean"
+
     if validate_input:
         validate_argument(locals(), "cluster_m", float, argmin=0, eqmin=True)
         validate_argument(locals(), "cluster_z", float, argmin=0, eqmin=True)
@@ -224,7 +228,7 @@ def generate_galaxy_catalog(
         validate_argument(locals(), "ngals", float, none_ok=True)
         validate_argument(locals(), "ngal_density", float, none_ok=True)
         validate_argument(locals(), "pz_bins", (int, "array"))
-        _validate_coordinate_system(locals(), "coordinate_system", str)
+        _validate_coordinate_system(locals(), "coordinate_system")
 
     if zsrc_min is None:
         zsrc_min = cluster_z + 0.1
@@ -393,7 +397,7 @@ def _generate_galaxy_catalog(
     # pylint: disable=R0914
 
     # Set the source galaxy redshifts
-    galaxy_catalog = _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals)
+    galaxy_catalog = _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals, coordinate_system)
 
     # Add photo-z errors and pdfs to source galaxy redshifts
     if photoz_sigma_unscaled is not None:
@@ -486,7 +490,7 @@ def _generate_galaxy_catalog(
     return galaxy_catalog[cols]
 
 
-def _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals):
+def _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals, coordinate_system):
     """Set source galaxy redshifts either set to a fixed value or draw from a predefined
     distribution. Return a table (GCData) of the source galaxies
 
@@ -542,7 +546,9 @@ def _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals):
     else:
         raise ValueError(f"zsrc must be a float, chang13 or desc_srd. You set: {zsrc}")
 
-    return GCData([zsrc_list, zsrc_list], names=("ztrue", "z"))
+    return GCData(
+        [zsrc_list, zsrc_list], names=("ztrue", "z"), meta={"coordinate_system": coordinate_system}
+    )
 
 
 def _compute_photoz_pdfs(

clmm/utils/validation.py

@@ -226,24 +226,21 @@ def _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c):
     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):
+
+def _validate_coordinate_system(loc, argname):
     r"""Validate the coordinate system.
 
     Parameters
     ----------
     loc: dict
         Dictionary with all input arguments. Should be locals().
-    coordinate_system: str
-        Coordinate system of the ellipticity components. Must be either 'celestial' or 'euclidean'.
-    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
+    argname: str
+        Name of argument to be tested.
     """
-    validate_argument(loc, coordinate_system, valid_type)
-    if loc[coordinate_system] not in ["celestial", "euclidean"]:
-        raise ValueError(f"{coordinate_system} must be 'celestial' or 'euclidean'.")
+    validate_argument(loc, argname, str)
+    if loc[argname] not in ["celestial", "euclidean"]:
+        raise ValueError(f"{argname} must be 'celestial' or 'euclidean'.")
+
 
 class DiffArray:
     """Array where arr1==arr2 is actually all(arr1==arr)"""

docs/doc-config.ini

@@ -33,6 +33,6 @@ OTHER
 ../examples/demo_compute_deltasigma_weights.ipynb
 ../examples/demo_mass_conversion.ipynb
 ../examples/demo_mock_ensemble_realistic.ipynb
-../examples/demo_coordinate_system_datasets.ipynb
+../examples/demo_coordinate_system.ipynb
 #../examples/other_compare_NFW_critical_massdef.ipynb
 #../examples/other_flat_sky_limitations.ipynb