Merge pull request yt-project#4839 from chrishavlin/geographic_spheri…

…cal_index_fixes

Bug: fix geographic coordinate conversions

nose_ignores.txt

@@ -12,6 +12,7 @@
 --ignore-file=test_file_sanitizer\.py
 --ignore-file=test_firefly\.py
 --ignore-file=test_geometries\.py
+--ignore-file=test_geographic_coordinates\.py
 --ignore-file=test_glue\.py
 --ignore-file=test_image_comp_2D_plots\.py
 --ignore-file=test_image_comp_geo\.py

tests/tests.yaml

@@ -183,6 +183,7 @@ other_tests:
      - "--ignore-file=test_file_sanitizer\\.py"
      - "--ignore-file=test_firefly\\.py"
      - "--ignore-file=test_geometries\\.py"
+     - "--ignore-file=test_geographic_coordinates\\.py"
      - "--ignore-file=test_glue\\.py"
      - "--ignore-file=test_image_comp_2D_plots\\.py"
      - "--ignore-file=test_image_comp_geo\\.py"

yt/geometry/coordinates/geographic_coordinates.py

@@ -125,7 +125,7 @@ def _path_longitude(field, data):
                 * data["index", "dlongitude"]
                 * np.pi
                 / 180.0
-                * np.sin((data["index", "latitude"] + 90.0) * np.pi / 180.0)
+                * np.sin((90 - data["index", "latitude"]) * np.pi / 180.0)
             )
 
         registry.add_field(
@@ -137,7 +137,8 @@ def _path_longitude(field, data):
 
         def _latitude_to_theta(field, data):
             # latitude runs from -90 to 90
-            return (data["index", "latitude"] + 90) * np.pi / 180.0
+            # theta = 0 at +90 deg, np.pi at -90
+            return (90.0 - data["index", "latitude"]) * np.pi / 180.0
 
         registry.add_field(
             ("index", "theta"),
@@ -158,7 +159,11 @@ def _dlatitude_to_dtheta(field, data):
 
         def _longitude_to_phi(field, data):
             # longitude runs from -180 to 180
-            return (data["index", "longitude"] + 180) * np.pi / 180.0
+            lonvals = data[("index", "longitude")]
+            neglons = lonvals < 0.0
+            if np.any(neglons):
+                lonvals[neglons] = lonvals[neglons] + 360.0
+            return lonvals * np.pi / 180.0
 
         registry.add_field(
             ("index", "phi"), sampling_type="cell", function=_longitude_to_phi, units=""
@@ -318,16 +323,16 @@ def convert_to_cartesian(self, coord):
             lon = self.axis_id["longitude"]
             lat = self.axis_id["latitude"]
             r = factor * coord[:, rad] + offset
-            theta = coord[:, lon] * np.pi / 180
-            phi = coord[:, lat] * np.pi / 180
+            theta = (90.0 - coord[:, lat]) * np.pi / 180
+            phi = coord[:, lon] * np.pi / 180
             nc = np.zeros_like(coord)
             # r, theta, phi
             nc[:, lat] = np.cos(phi) * np.sin(theta) * r
             nc[:, lon] = np.sin(phi) * np.sin(theta) * r
             nc[:, rad] = np.cos(theta) * r
         else:
             a, b, c = coord
-            theta = b * np.pi / 180
+            theta = (90.0 - b) * np.pi / 180
             phi = a * np.pi / 180
             r = factor * c + offset
             nc = (

yt/geometry/coordinates/tests/test_geographic_coordinates.py

@@ -1,6 +1,7 @@
 # Some tests for the geographic coordinates handler
 
 import numpy as np
+import pytest
 from numpy.testing import assert_equal
 
 from yt.testing import assert_rel_equal, fake_amr_ds
@@ -9,16 +10,25 @@
 # compute our volume correctly.
 
 
-def test_geographic_coordinates():
+@pytest.mark.parametrize("geometry", ("geographic", "internal_geographic"))
+def test_geographic_coordinates(geometry):
     # We're going to load up a simple AMR grid and check its volume
     # calculations and path length calculations.
 
     # Note that we are setting it up to have an altitude of 1000 maximum, which
     # means our volume will be that of a shell 1000 wide, starting at r of
     # whatever our surface_height is set to.
-    ds = fake_amr_ds(geometry="geographic")
-    ds.surface_height = ds.quan(5000.0, "code_length")
-    axes = ["latitude", "longitude", "altitude"]
+    ds = fake_amr_ds(geometry=geometry)
+    if geometry == "geographic":
+        ds.surface_height = ds.quan(5000.0, "code_length")
+        inner_r = ds.surface_height
+        outer_r = ds.surface_height + ds.domain_width[2]
+    else:
+        ds.outer_radius = ds.quan(5000.0, "code_length")
+        inner_r = ds.outer_radius - ds.domain_right_edge[2]
+        outer_r = ds.outer_radius
+    radial_axis = ds.coordinates.radial_axis
+    axes = ds.coordinates.axis_order
     for i, axis in enumerate(axes):
         dd = ds.all_data()
         fi = ("index", axis)
@@ -28,18 +38,19 @@ def test_geographic_coordinates():
         mi = np.argmin(dd[fi])
         assert_equal(dd[fi][mi] - dd[fd][mi] / 2.0, ds.domain_left_edge[i].d)
         assert_equal(dd[fd].max(), (ds.domain_width / ds.domain_dimensions)[i].d)
-    inner_r = ds.surface_height
-    outer_r = ds.surface_height + ds.domain_width[2]
+
     assert_equal(dd["index", "dtheta"], dd["index", "dlatitude"] * np.pi / 180.0)
     assert_equal(dd["index", "dphi"], dd["index", "dlongitude"] * np.pi / 180.0)
     # Note our terrible agreement here.
     assert_rel_equal(
         dd["index", "cell_volume"].sum(dtype="float64"),
         (4.0 / 3.0) * np.pi * (outer_r**3 - inner_r**3),
-        10,
+        14,
+    )
+    assert_equal(
+        dd["index", f"path_element_{radial_axis}"], dd["index", f"d{radial_axis}"]
     )
-    assert_equal(dd["index", "path_element_altitude"], dd["index", "daltitude"])
-    assert_equal(dd["index", "path_element_altitude"], dd["index", "dr"])
+    assert_equal(dd["index", f"path_element_{radial_axis}"], dd["index", "dr"])
     # Note that latitude corresponds to theta, longitude to phi
     assert_equal(
         dd["index", "path_element_latitude"],
@@ -52,57 +63,58 @@ def test_geographic_coordinates():
             * dd["index", "dlongitude"]
             * np.pi
             / 180.0
-            * np.sin((dd["index", "latitude"] + 90.0) * np.pi / 180.0)
+            * np.sin((90 - dd["index", "latitude"]) * np.pi / 180.0)
         ),
     )
     # We also want to check that our radius is correct
-    assert_equal(dd["index", "r"], dd["index", "altitude"] + ds.surface_height)
+    offset, factor = ds.coordinates._retrieve_radial_offset()
+    radius = factor * dd["index", radial_axis] + offset
+    assert_equal(dd["index", "r"], radius)
 
 
-def test_internal_geographic_coordinates():
-    # We're going to load up a simple AMR grid and check its volume
-    # calculations and path length calculations.
+@pytest.mark.parametrize("geometry", ("geographic", "internal_geographic"))
+def test_geographic_conversions(geometry):
 
-    # Note that we are setting it up to have depth of 1000 maximum, which
-    # means our volume will be that of a shell 1000 wide, starting at r of
-    # outer_radius - 1000.
-    ds = fake_amr_ds(geometry="internal_geographic")
-    ds.outer_radius = ds.quan(5000.0, "code_length")
-    axes = ["latitude", "longitude", "depth"]
-    for i, axis in enumerate(axes):
-        dd = ds.all_data()
-        fi = ("index", axis)
-        fd = ("index", f"d{axis}")
-        ma = np.argmax(dd[fi])
-        assert_equal(dd[fi][ma] + dd[fd][ma] / 2.0, ds.domain_right_edge[i].d)
-        mi = np.argmin(dd[fi])
-        assert_equal(dd[fi][mi] - dd[fd][mi] / 2.0, ds.domain_left_edge[i].d)
-        assert_equal(dd[fd].max(), (ds.domain_width / ds.domain_dimensions)[i].d)
-    inner_r = ds.outer_radius - ds.domain_right_edge[2]
-    outer_r = ds.outer_radius
-    assert_equal(dd["index", "dtheta"], dd["index", "dlatitude"] * np.pi / 180.0)
-    assert_equal(dd["index", "dphi"], dd["index", "dlongitude"] * np.pi / 180.0)
-    assert_rel_equal(
-        dd["index", "cell_volume"].sum(dtype="float64"),
-        (4.0 / 3.0) * np.pi * (outer_r**3 - inner_r**3),
-        10,
-    )
-    assert_equal(dd["index", "path_element_depth"], dd["index", "ddepth"])
-    assert_equal(dd["index", "path_element_depth"], dd["index", "dr"])
-    # Note that latitude corresponds to theta, longitude to phi
-    assert_equal(
-        dd["index", "path_element_latitude"],
-        dd["index", "r"] * dd["index", "dlatitude"] * np.pi / 180.0,
-    )
-    assert_equal(
-        dd["index", "path_element_longitude"],
-        (
-            dd["index", "r"]
-            * dd["index", "dlongitude"]
-            * np.pi
-            / 180.0
-            * np.sin((dd["index", "latitude"] + 90.0) * np.pi / 180.0)
-        ),
-    )
-    # We also want to check that our radius is correct
-    assert_equal(dd["index", "r"], -1.0 * dd["index", "depth"] + ds.outer_radius)
+    ds = fake_amr_ds(geometry=geometry)
+    ad = ds.all_data()
+    lats = ad["index", "latitude"]
+    dlats = ad["index", "dlatitude"]
+    theta = ad["index", "theta"]
+    dtheta = ad["index", "dtheta"]
+
+    # check that theta = 0, pi at latitudes of 90, -90
+    south_pole_id = np.where(lats == np.min(lats))[0][0]
+    north_pole_id = np.where(lats == np.max(lats))[0][0]
+
+    # check that we do in fact have -90, 90 exactly
+    assert lats[south_pole_id] - dlats[south_pole_id] / 2.0 == -90.0
+    assert lats[north_pole_id] + dlats[north_pole_id] / 2.0 == 90.0
+
+    # check that theta=0 at the north pole, np.pi at the south
+    assert theta[north_pole_id] - dtheta[north_pole_id] / 2 == 0.0
+    assert theta[south_pole_id] + dtheta[south_pole_id] / 2 == np.pi
+
+    # check that longitude-phi conversions
+    phi = ad["index", "phi"]
+    dphi = ad["index", "dphi"]
+    lons = ad["index", "longitude"]
+    dlon = ad["index", "dlongitude"]
+    lon_180 = np.where(lons == np.max(lons))[0][0]
+    lon_neg180 = np.where(lons == np.min(lons))[0][0]
+    # check we have -180, 180 exactly
+    assert lons[lon_neg180] - dlon[lon_neg180] / 2.0 == -180.0
+    assert lons[lon_180] + dlon[lon_180] / 2.0 == 180.0
+    # check that those both convert to phi = np.pi
+    assert phi[lon_neg180] - dphi[lon_neg180] / 2.0 == np.pi
+    assert phi[lon_180] + dphi[lon_180] / 2.0 == np.pi
+
+    # check that z = +/- radius at +/-90
+    # default expected axis order: lat, lon, radial axis
+    r_val = ds.coordinates._retrieve_radial_offset()[0]
+    coords = np.zeros((2, 3))
+    coords[0, 0] = 90.0
+    coords[1, 0] = -90.0
+    xyz = ds.coordinates.convert_to_cartesian(coords)
+    z = xyz[:, 2]
+    assert z[0] == r_val
+    assert z[1] == -r_val