diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py
index 860cee6c2..fff56070a 100644
--- a/clmm/dataops/__init__.py
+++ b/clmm/dataops/__init__.py
@@ -378,13 +378,13 @@ def _compute_lensing_angles_flatsky(ra_lens, dec_lens, ra_source_list, dec_sourc
     # deltax[deltax < -np.pi] = deltax[deltax < -np.pi]+2.*np.pi
     angsep = np.sqrt(deltax**2 + deltay**2)
     phi = np.arctan2(deltay, -deltax)
-    if coordinate_system == "sky":
-        phi = np.pi - phi
     # Forcing phi to be zero everytime angsep is zero. This is necessary due to arctan2 features.
     if np.iterable(phi):
         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
@@ -417,14 +417,14 @@ def _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_list, dec_sourc
     angsep, phi = sk_lens.separation(sk_src).rad, sk_lens.position_angle(sk_src).rad
     # Transformations for phi to have same orientation as _compute_lensing_angles_flatsky
     phi += 0.5 * np.pi
-    if coordinate_system == "sky":
-        phi = np.pi - phi
     if np.iterable(phi):
         phi[phi > np.pi] -= 2 * np.pi
         phi[angsep == 0] = 0
     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
 
 
diff --git a/tests/test_dataops.py b/tests/test_dataops.py
index 9514a3905..192bba03b 100644
--- a/tests/test_dataops.py
+++ b/tests/test_dataops.py
@@ -155,6 +155,31 @@ def test_compute_lensing_angles_flatsky():
     )
 
 
+def test_compute_lensing_angles_astropy():
+    """test compute lensing angles astropy"""
+
+    # coordinate_system conversion
+    ra_l, dec_l = 161.32, 51.49
+    ra_s, dec_s = np.array([161.29, 161.34]), np.array([51.45, 51.55])
+    thetas_pixel, phis_pixel = da._compute_lensing_angles_astropy(ra_l, dec_l, ra_s, dec_s, coordinate_system="pixel")
+    thetas_sky, phis_sky = da._compute_lensing_angles_astropy(ra_l, dec_l, ra_s, dec_s, coordinate_system="sky")
+
+    print(phis_pixel, phis_sky)
+
+    assert_allclose(
+        thetas_sky,
+        thetas_pixel,
+        **TOLERANCE,
+        err_msg="Conversion from sky to pixel coordinate system for theta failed",
+    )
+
+    assert_allclose(
+        phis_sky,
+        np.pi - phis_pixel,
+        **TOLERANCE,
+        err_msg="Conversion from sky to pixel coordinate system for phi failed",
+    )
+
 
 def test_compute_tangential_and_cross_components(modeling_data):
     """test compute tangential and cross components"""