Catalog crossmatch example runs into an error with units
Running the exact example code provided in https://lscsoft.docs.ligo.org/ligo.skymap/postprocess/crossmatch.html, I run into this error:
UnitConversionError: 'Unit("Mpc2")' is not a scaled version of 'Unit("Mpc")'
and
UnitConversionError: 'Mpc2' (area) and 'Mpc' (length) are not convertible
.
Here is the full traceback:
10000it [00:00, 41431.73it/s]
---------------------------------------------------------------------------
UnitConversionError Traceback (most recent call last)
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/quantity.py:982, in Quantity.to_value(self, unit, equivalencies)
981 try:
--> 982 scale = self.unit._to(unit)
983 except Exception:
984 # Short-cut failed; try default (maybe equivalencies help).
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/core.py:1159, in UnitBase._to(self, other)
1157 return self_decomposed.scale / other_decomposed.scale
-> 1159 raise UnitConversionError(f"'{self!r}' is not a scaled version of '{other!r}'")
UnitConversionError: 'Unit("Mpc2")' is not a scaled version of 'Unit("Mpc")'
During handling of the above exception, another exception occurred:
UnitConversionError Traceback (most recent call last)
Cell In[82], line 18
15 url = 'https://gracedb.ligo.org/api/superevents/S190814bv/files/bayestar.multiorder.fits'
16 skymap = read_sky_map(url, moc=True)
---> 18 result = crossmatch(ex_skymap, gal_coords)
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/ligo/skymap/postprocess/crossmatch.py:295, in crossmatch(sky_map, coordinates, contours, areas, modes, cosmology)
293 true_dec = coordinates.dec.rad
294 if np.any(coordinates.distance != 1):
--> 295 true_dist = coordinates.distance.to_value(u.Mpc)
296 else:
297 true_dist = None
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/quantity.py:985, in Quantity.to_value(self, unit, equivalencies)
982 scale = self.unit._to(unit)
983 except Exception:
984 # Short-cut failed; try default (maybe equivalencies help).
--> 985 value = self._to_value(unit, equivalencies)
986 else:
987 value = self.view(np.ndarray)
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/quantity.py:891, in Quantity._to_value(self, unit, equivalencies)
888 equivalencies = self._equivalencies
889 if not self.dtype.names or isinstance(self.unit, StructuredUnit):
890 # Standard path, let unit to do work.
--> 891 return self.unit.to(
892 unit, self.view(np.ndarray), equivalencies=equivalencies
893 )
895 else:
896 # The .to() method of a simple unit cannot convert a structured
897 # dtype, so we work around it, by recursing.
898 # TODO: deprecate this?
899 # Convert simple to Structured on initialization?
900 result = np.empty_like(self.view(np.ndarray))
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/core.py:1195, in UnitBase.to(self, other, value, equivalencies)
1193 return UNITY
1194 else:
-> 1195 return self._get_converter(Unit(other), equivalencies)(value)
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/core.py:1124, in UnitBase._get_converter(self, other, equivalencies)
1121 else:
1122 return lambda v: b(converter(v))
-> 1124 raise exc
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/core.py:1107, in UnitBase._get_converter(self, other, equivalencies)
1105 # if that doesn't work, maybe we can do it with equivalencies?
1106 try:
-> 1107 return self._apply_equivalencies(
1108 self, other, self._normalize_equivalencies(equivalencies)
1109 )
1110 except UnitsError as exc:
1111 # Last hope: maybe other knows how to do it?
1112 # We assume the equivalencies have the unit itself as first item.
1113 # TODO: maybe better for other to have a `_back_converter` method?
1114 if hasattr(other, "equivalencies"):
File ~/miniconda3/envs/catalog-pe/lib/python3.11/site-packages/astropy/units/core.py:1085, in UnitBase._apply_equivalencies(self, unit, other, equivalencies)
1082 unit_str = get_err_str(unit)
1083 other_str = get_err_str(other)
-> 1085 raise UnitConversionError(f"{unit_str} and {other_str} are not convertible")
UnitConversionError: 'Mpc2' (area) and 'Mpc' (length) are not convertible