From ea9bbbeda7494ba67b1d8d07910599a6f9ffca11 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Tue, 5 Sep 2023 15:16:13 +1000
Subject: [PATCH 01/23] simplify angle_between

---
 regional_mom6/regional_mom6.py | 25 +++++++------------------
 tests/test_grid_generation.py  |  7 +++++++
 2 files changed, 14 insertions(+), 18 deletions(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index 0818c14a..f5719a24 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -261,7 +261,6 @@ def dz(npoints, ratio, target_depth, min_dz=0.0001, tolerance=1):
 
     Returns:
         numpy.array: An array containing the thickness profile.
-
     """
 
     profile = min_dz + 0.5 * (np.abs(ratio) * min_dz - min_dz) * (
@@ -279,24 +278,14 @@ def dz(npoints, ratio, target_depth, min_dz=0.0001, tolerance=1):
     return dz(npoints, ratio, target_depth, min_dz * err_ratio)
 
 
-# Borrowed from grid tools (GFDL)
 def angle_between(v1, v2, v3):
-    """Returns angle v2-v1-v3 i.e betweeen v1-v2 and v1-v3."""
-
-    # vector product between v1 and v2
-    px = v1[1] * v2[2] - v1[2] * v2[1]
-    py = v1[2] * v2[0] - v1[0] * v2[2]
-    pz = v1[0] * v2[1] - v1[1] * v2[0]
-    # vector product between v1 and v3
-    qx = v1[1] * v3[2] - v1[2] * v3[1]
-    qy = v1[2] * v3[0] - v1[0] * v3[2]
-    qz = v1[0] * v3[1] - v1[1] * v3[0]
-
-    ddd = (px * px + py * py + pz * pz) * (qx * qx + qy * qy + qz * qz)
-    ddd = (px * qx + py * qy + pz * qz) / np.sqrt(ddd)
-    angle = np.arccos(ddd)
-
-    return angle
+    """Returns angle v2-v1-v3 that is between the vectors v1-v2 and v1-v3."""
+
+    v1xv2 = np.cross(v1, v2)
+    v1xv3 = np.cross(v1, v3)
+    cosangle = np.dot(v1xv2, v1xv3) / np.sqrt(np.dot(v1xv2, v1xv2) * np.dot(v1xv3, v1xv3))
+
+    return np.arccos(cosangle)
 
 
 # Borrowed from grid tools (GFDL)
diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index d39bed9e..b0b63a9e 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -5,6 +5,13 @@
 from regional_mom6 import rectangular_hgrid
 import xarray as xr
 
+def random_unit_vector():
+    """Return a random unit vector on the unit sphere."""
+    z = 2 * np.random.rand(1)[0] - 1; z
+    θ = 2 * np.pi * np.random.rand(1)[0]
+    return [np.sqrt(1 - z**2) * np.cos(θ), np.sqrt(1 - z**2) * np.sin(θ), z]
+
+
 
 @pytest.mark.parametrize(
     ("v1", "v2", "v3", "true_angle"),

From 812a302f20be965ac7860791d263329cb5134ff7 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Tue, 5 Sep 2023 15:21:39 +1000
Subject: [PATCH 02/23] black format

---
 regional_mom6/regional_mom6.py | 4 +++-
 tests/test_grid_generation.py  | 5 +++--
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index f5719a24..cf657876 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -283,7 +283,9 @@ def angle_between(v1, v2, v3):
 
     v1xv2 = np.cross(v1, v2)
     v1xv3 = np.cross(v1, v3)
-    cosangle = np.dot(v1xv2, v1xv3) / np.sqrt(np.dot(v1xv2, v1xv2) * np.dot(v1xv3, v1xv3))
+    cosangle = np.dot(v1xv2, v1xv3) / np.sqrt(
+        np.dot(v1xv2, v1xv2) * np.dot(v1xv3, v1xv3)
+    )
 
     return np.arccos(cosangle)
 
diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index b0b63a9e..4c50148e 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -5,14 +5,15 @@
 from regional_mom6 import rectangular_hgrid
 import xarray as xr
 
+
 def random_unit_vector():
     """Return a random unit vector on the unit sphere."""
-    z = 2 * np.random.rand(1)[0] - 1; z
+    z = 2 * np.random.rand(1)[0] - 1
+    z
     θ = 2 * np.pi * np.random.rand(1)[0]
     return [np.sqrt(1 - z**2) * np.cos(θ), np.sqrt(1 - z**2) * np.sin(θ), z]
 
 
-
 @pytest.mark.parametrize(
     ("v1", "v2", "v3", "true_angle"),
     [

From 3d3179ac21c36934be3e3ee4bf68fff22211e100 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Tue, 5 Sep 2023 16:04:40 +1000
Subject: [PATCH 03/23] remove random unit vector; add 2 more tests

---
 tests/test_grid_generation.py | 10 ++--------
 1 file changed, 2 insertions(+), 8 deletions(-)

diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 4c50148e..405d2c27 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -6,19 +6,13 @@
 import xarray as xr
 
 
-def random_unit_vector():
-    """Return a random unit vector on the unit sphere."""
-    z = 2 * np.random.rand(1)[0] - 1
-    z
-    θ = 2 * np.pi * np.random.rand(1)[0]
-    return [np.sqrt(1 - z**2) * np.cos(θ), np.sqrt(1 - z**2) * np.sin(θ), z]
-
-
 @pytest.mark.parametrize(
     ("v1", "v2", "v3", "true_angle"),
     [
         ([1, 0, 0], [0, 1, 0], [0, 0, 1], np.pi / 2),
         ([1, 0, 0], [1, 1, 0], [0, 1, 1], np.pi / 4),
+        ([1, 0, 0], [1, 1, 1], [0, 0, 1], np.pi / 4),
+        ([1, 1, 1], [1, 1, 0], [0, 1, 1], 2 * np.pi / 3)
     ],
 )
 def test_angle_between(v1, v2, v3, true_angle):

From 6721c471f80127079604901bb73bf69910e60d4f Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 08:35:41 +1000
Subject: [PATCH 04/23] refactor quadrilateral_area

---
 regional_mom6/regional_mom6.py | 39 +++++++++++++++++++++++++++-------
 tests/test_grid_generation.py  |  5 +++--
 2 files changed, 34 insertions(+), 10 deletions(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index cf657876..1c752a29 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -22,6 +22,7 @@
     "dz",
     "angle_between",
     "quadilateral_area",
+    "quadilateral_areas",
     "rectangular_hgrid",
     "experiment",
     "segment",
@@ -290,9 +291,23 @@ def angle_between(v1, v2, v3):
     return np.arccos(cosangle)
 
 
+def quadilateral_area(v1, v2, v3, v4):
+    """Returns area of a spherical quadrilateral on the unit sphere that
+    has vertices on 3-vectors `v1`, `v2`, `v3`, `v4` (counter-clockwise
+    orientation is implied). The area is computed as the excess of the
+    sum of the spherical angles of the quadrilateral from 2π."""
+
+    a1 = angle_between(v1, v2, v4)
+    a2 = angle_between(v2, v3, v1)
+    a3 = angle_between(v3, v4, v2)
+    a4 = angle_between(v4, v1, v3)
+
+    return a1 + a2 + a3 + a4 - 2.0 * np.pi
+
 # Borrowed from grid tools (GFDL)
-def quadilateral_area(lat, lon):
-    """Returns area of spherical quadrilaterals (bounded by great arcs)."""
+def quadilateral_areas(lat, lon):
+    """Returns area of spherical quadrilaterals on the unit sphere that are formed
+    by constant latitude and longitude lines on the `lat`-`lon` grid provided."""
 
     # x, y, z are 3D coordinates on the unit sphere
     x = np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(lon))
@@ -304,12 +319,20 @@ def quadilateral_area(lat, lon):
     c2 = (x[1:, 1:], y[1:, 1:], z[1:, 1:])
     c3 = (x[1:, :-1], y[1:, :-1], z[1:, :-1])
 
-    a0 = angle_between(c1, c0, c2)
-    a1 = angle_between(c2, c1, c3)
-    a2 = angle_between(c3, c2, c0)
-    a3 = angle_between(c0, c3, c1)
+    nx, ny = np.shape(lat)
+
+    areas = np.zeros((nx-1, ny-1))
+
+    for j in range(ny-1):
+        for i in range(nx-1):
+            v1 = [x[ i,   j ], y[ i,   j ], z[ i,   j ]]
+            v2 = [x[ i,  j+1], y[ i,  j+1], z[ i,  j+1]]
+            v3 = [x[i+1, j+1], y[i+1, j+1], z[i+1, j+1]]
+            v4 = [x[i+1,  j ], y[i+1,  j ], z[i+1,  j ]]
+
+            areas[i, j] = quadilateral_area(v1, v2, v3, v4)
 
-    return a0 + a1 + a2 + a3 - 2.0 * np.pi
+    return areas
 
 
 def rectangular_hgrid(λ, φ):
@@ -347,7 +370,7 @@ def rectangular_hgrid(λ, φ):
 
     lon, lat = np.meshgrid(λ, φ)
 
-    area = quadilateral_area(lat, lon) * R**2
+    area = quadilateral_areas(lat, lon) * R**2
 
     attrs = {
         "tile": {
diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 405d2c27..4cfcfe9b 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -2,6 +2,7 @@
 import pytest
 from regional_mom6 import angle_between
 from regional_mom6 import quadilateral_area
+from regional_mom6 import quadilateral_areas
 from regional_mom6 import rectangular_hgrid
 import xarray as xr
 
@@ -33,8 +34,8 @@ def test_angle_between(v1, v2, v3, true_angle):
         (lat2, lon2, np.pi),
     ],
 )
-def test_quadilateral_area(lat, lon, true_area):
-    assert np.isclose(np.sum(quadilateral_area(lat, lon)), true_area)
+def test_quadilateral_areas(lat, lon, true_area):
+    assert np.isclose(np.sum(quadilateral_areas(lat, lon)), true_area)
 
 
 # a simple test that rectangular_hgrid runs without erroring

From d414ea17e36425bb87a466981db030fa2f7d4c0b Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 08:36:03 +1000
Subject: [PATCH 05/23] black format

---
 regional_mom6/regional_mom6.py | 15 ++++++++-------
 tests/test_grid_generation.py  |  2 +-
 2 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index 1c752a29..1b1b9100 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -304,6 +304,7 @@ def quadilateral_area(v1, v2, v3, v4):
 
     return a1 + a2 + a3 + a4 - 2.0 * np.pi
 
+
 # Borrowed from grid tools (GFDL)
 def quadilateral_areas(lat, lon):
     """Returns area of spherical quadrilaterals on the unit sphere that are formed
@@ -321,14 +322,14 @@ def quadilateral_areas(lat, lon):
 
     nx, ny = np.shape(lat)
 
-    areas = np.zeros((nx-1, ny-1))
+    areas = np.zeros((nx - 1, ny - 1))
 
-    for j in range(ny-1):
-        for i in range(nx-1):
-            v1 = [x[ i,   j ], y[ i,   j ], z[ i,   j ]]
-            v2 = [x[ i,  j+1], y[ i,  j+1], z[ i,  j+1]]
-            v3 = [x[i+1, j+1], y[i+1, j+1], z[i+1, j+1]]
-            v4 = [x[i+1,  j ], y[i+1,  j ], z[i+1,  j ]]
+    for j in range(ny - 1):
+        for i in range(nx - 1):
+            v1 = [x[i, j], y[i, j], z[i, j]]
+            v2 = [x[i, j + 1], y[i, j + 1], z[i, j + 1]]
+            v3 = [x[i + 1, j + 1], y[i + 1, j + 1], z[i + 1, j + 1]]
+            v4 = [x[i + 1, j], y[i + 1, j], z[i + 1, j]]
 
             areas[i, j] = quadilateral_area(v1, v2, v3, v4)
 
diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 4cfcfe9b..7b37ba78 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -13,7 +13,7 @@
         ([1, 0, 0], [0, 1, 0], [0, 0, 1], np.pi / 2),
         ([1, 0, 0], [1, 1, 0], [0, 1, 1], np.pi / 4),
         ([1, 0, 0], [1, 1, 1], [0, 0, 1], np.pi / 4),
-        ([1, 1, 1], [1, 1, 0], [0, 1, 1], 2 * np.pi / 3)
+        ([1, 1, 1], [1, 1, 0], [0, 1, 1], 2 * np.pi / 3),
     ],
 )
 def test_angle_between(v1, v2, v3, true_angle):

From 9ad060ed3f6d01ca1d8bb411d871c87773ccaf8c Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 08:41:33 +1000
Subject: [PATCH 06/23] adds comment

---
 regional_mom6/regional_mom6.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index bc82834a..51ed3db3 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -321,6 +321,7 @@ def quadilateral_areas(lat, lon):
 
     for j in range(ny - 1):
         for i in range(nx - 1):
+            # construct the 4 3-vectors v1, v2, v3, v4 that point to the vertices of the quadrilaterals
             v1 = [x[i, j], y[i, j], z[i, j]]
             v2 = [x[i, j + 1], y[i, j + 1], z[i, j + 1]]
             v3 = [x[i + 1, j + 1], y[i + 1, j + 1], z[i + 1, j + 1]]

From ad69968917089b738ed6ea791a3df1466370a647 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 08:47:49 +1000
Subject: [PATCH 07/23] better docstring for quadilateral_areas

---
 regional_mom6/regional_mom6.py | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index 51ed3db3..22fd2628 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -305,10 +305,19 @@ def quadilateral_area(v1, v2, v3, v4):
     return a1 + a2 + a3 + a4 - 2.0 * np.pi
 
 
-# Borrowed from grid tools (GFDL)
 def quadilateral_areas(lat, lon):
     """Returns area of spherical quadrilaterals on the unit sphere that are formed
-    by constant latitude and longitude lines on the `lat`-`lon` grid provided."""
+    by constant latitude and longitude lines on the `lat`-`lon` grid provided.
+
+    Args:
+        lat (array): Array of latitude points
+        lon (array): Array of longitude points
+
+    Returns:
+        areas (array): Array with the areas of the quadrilaterals defined by the
+                       `lat`-`lon` grid provided. If the `lat`-`lon` are `m x n`
+                       then `areas` is `(m-1) x (n-1)`.
+    """
 
     # x, y, z are 3D coordinates on the unit sphere
     x = np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(lon))

From d34d2f8112103c8a2a285f0b6ae8a3c2e885cfbc Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 08:58:19 +1000
Subject: [PATCH 08/23] add test for quadilateral_area

---
 tests/test_grid_generation.py | 26 ++++++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 7b37ba78..5df9c921 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -20,6 +20,32 @@ def test_angle_between(v1, v2, v3, true_angle):
     assert np.isclose(angle_between(v1, v2, v3), true_angle)
 
 
+def latlon_to_cartesian(lat, lon):
+    """Convert latitude-longitude to Cartesian coordinates on a unit sphere."""
+    x = np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(lon))
+    y = np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(lon))
+    z = np.sin(np.deg2rad(lat))
+    return np.array([x, y, z])
+
+
+@pytest.mark.parametrize(
+    ("v1", "v2", "v3", "v4", "true_area"),
+    [
+        (
+            latlon_to_cartesian(0, 0),
+            latlon_to_cartesian(0, 90),
+            latlon_to_cartesian(90, 0),
+            latlon_to_cartesian(0, -90),
+            np.pi,
+        ),
+    ],
+)
+def test_quadilateral_area(v1, v2, v3, v4, true_area):
+    print(quadilateral_area(v1, v2, v3, v4))
+    print(true_area)
+    assert np.isclose(quadilateral_area(v1, v2, v3, v4), true_area)
+
+
 # create a lat-lon mesh that covers 1/4 of the North Hemisphere
 lon1, lat1 = np.meshgrid(np.linspace(0, 90, 5), np.linspace(0, 90, 5))
 

From 283778a169e6b7e6e8cb9486d7b93083c08ba3ec Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 09:21:03 +1000
Subject: [PATCH 09/23] add latlon_to_cartesian + test

---
 regional_mom6/regional_mom6.py | 14 ++++++++++----
 tests/test_grid_generation.py  | 22 ++++++++++++++--------
 2 files changed, 24 insertions(+), 12 deletions(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index 22fd2628..c8564f07 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -21,6 +21,7 @@
     "motu_requests",
     "dz",
     "angle_between",
+    "latlon_to_cartesian",
     "quadilateral_area",
     "quadilateral_areas",
     "rectangular_hgrid",
@@ -305,6 +306,14 @@ def quadilateral_area(v1, v2, v3, v4):
     return a1 + a2 + a3 + a4 - 2.0 * np.pi
 
 
+def latlon_to_cartesian(lat, lon):
+    """Convert latitude-longitude to Cartesian coordinates on a unit sphere."""
+    x = np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(lon))
+    y = np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(lon))
+    z = np.sin(np.deg2rad(lat))
+    return x, y, z
+
+
 def quadilateral_areas(lat, lon):
     """Returns area of spherical quadrilaterals on the unit sphere that are formed
     by constant latitude and longitude lines on the `lat`-`lon` grid provided.
@@ -319,10 +328,7 @@ def quadilateral_areas(lat, lon):
                        then `areas` is `(m-1) x (n-1)`.
     """
 
-    # x, y, z are 3D coordinates on the unit sphere
-    x = np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(lon))
-    y = np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(lon))
-    z = np.sin(np.deg2rad(lat))
+    x, y, z = latlon_to_cartesian(lat, lon)
 
     nx, ny = np.shape(lat)
 
diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 5df9c921..3d1b1952 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -1,12 +1,26 @@
 import numpy as np
 import pytest
 from regional_mom6 import angle_between
+from regional_mom6 import latlon_to_cartesian
 from regional_mom6 import quadilateral_area
 from regional_mom6 import quadilateral_areas
 from regional_mom6 import rectangular_hgrid
 import xarray as xr
 
 
+@pytest.mark.parametrize(
+    ("lat", "lon", "true_xyz"),
+    [
+        (0, 0, (1, 0, 0)),
+        (90, 0, (0, 0, 1)),
+        (0, 90, (0, 1, 0)),
+        (-90, 0, (0, 0, -1)),
+    ],
+)
+def test_latlon_to_cartesian(lat, lon, true_xyz):
+    assert np.isclose(latlon_to_cartesian(lat, lon), true_xyz).all()
+
+
 @pytest.mark.parametrize(
     ("v1", "v2", "v3", "true_angle"),
     [
@@ -20,14 +34,6 @@ def test_angle_between(v1, v2, v3, true_angle):
     assert np.isclose(angle_between(v1, v2, v3), true_angle)
 
 
-def latlon_to_cartesian(lat, lon):
-    """Convert latitude-longitude to Cartesian coordinates on a unit sphere."""
-    x = np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(lon))
-    y = np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(lon))
-    z = np.sin(np.deg2rad(lat))
-    return np.array([x, y, z])
-
-
 @pytest.mark.parametrize(
     ("v1", "v2", "v3", "v4", "true_area"),
     [

From 4cf6a3b2626b4500c2ed47411d58cb19442e4055 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 09:23:53 +1000
Subject: [PATCH 10/23] one more test for quadilateral_area

---
 tests/test_grid_generation.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 3d1b1952..77fe2af1 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -44,6 +44,13 @@ def test_angle_between(v1, v2, v3, true_angle):
             latlon_to_cartesian(0, -90),
             np.pi,
         ),
+        (
+            latlon_to_cartesian(0, 0),
+            latlon_to_cartesian(90, 0),
+            latlon_to_cartesian(0, 90),
+            latlon_to_cartesian(-90, 0),
+            np.pi,
+        ),
     ],
 )
 def test_quadilateral_area(v1, v2, v3, v4, true_area):

From d4f8668c5a8da18fff701ffb4af2598ec1dbc0a9 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 09:24:25 +1000
Subject: [PATCH 11/23] remove debuggint print statements

---
 tests/test_grid_generation.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 77fe2af1..26375b44 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -54,8 +54,6 @@ def test_angle_between(v1, v2, v3, true_angle):
     ],
 )
 def test_quadilateral_area(v1, v2, v3, v4, true_area):
-    print(quadilateral_area(v1, v2, v3, v4))
-    print(true_area)
     assert np.isclose(quadilateral_area(v1, v2, v3, v4), true_area)
 
 

From 2366869774f73e5086a488cccff0b2d036f70000 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 09:25:55 +1000
Subject: [PATCH 12/23] cleaner test_quadilateral_areas

---
 tests/test_grid_generation.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 26375b44..6dda9066 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -59,16 +59,18 @@ def test_quadilateral_area(v1, v2, v3, v4, true_area):
 
 # create a lat-lon mesh that covers 1/4 of the North Hemisphere
 lon1, lat1 = np.meshgrid(np.linspace(0, 90, 5), np.linspace(0, 90, 5))
+area1 = 1 / 8 * (4 * np.pi)
 
 # create a lat-lon mesh that covers 1/4 of the whole globe
 lon2, lat2 = np.meshgrid(np.linspace(-45, 45, 5), np.linspace(-90, 90, 5))
+area2 = 1 / 4 * (4 * np.pi)
 
 
 @pytest.mark.parametrize(
     ("lat", "lon", "true_area"),
     [
-        (lat1, lon1, 0.5 * np.pi),
-        (lat2, lon2, np.pi),
+        (lat1, lon1, area1),
+        (lat2, lon2, area2),
     ],
 )
 def test_quadilateral_areas(lat, lon, true_area):

From b5877f33638687087571b4e3c6dff0418d0521ab Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 09:32:01 +1000
Subject: [PATCH 13/23] add clarification in latlon_to_cartesian doc

---
 regional_mom6/regional_mom6.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index c8564f07..a44120fa 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -307,7 +307,7 @@ def quadilateral_area(v1, v2, v3, v4):
 
 
 def latlon_to_cartesian(lat, lon):
-    """Convert latitude-longitude to Cartesian coordinates on a unit sphere."""
+    """Convert latitude-longitude (in degrees) to Cartesian coordinates on a unit sphere."""
     x = np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(lon))
     y = np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(lon))
     z = np.sin(np.deg2rad(lat))

From 00ad0b58f7683ec6d3b749541ba7c940d1c76083 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 09:32:19 +1000
Subject: [PATCH 14/23] add clarification in latlon_to_cartesian doc

---
 regional_mom6/regional_mom6.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index a44120fa..e3547130 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -319,8 +319,8 @@ def quadilateral_areas(lat, lon):
     by constant latitude and longitude lines on the `lat`-`lon` grid provided.
 
     Args:
-        lat (array): Array of latitude points
-        lon (array): Array of longitude points
+        lat (array): Array of latitude points (in degrees)
+        lon (array): Array of longitude points (in degrees)
 
     Returns:
         areas (array): Array with the areas of the quadrilaterals defined by the

From 2ee793d88c9a817419f2f5bd45fb6b34c9452b8c Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 09:46:39 +1000
Subject: [PATCH 15/23] add clarification in angle_between docs

---
 regional_mom6/regional_mom6.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index e3547130..43a626d4 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -281,10 +281,11 @@ def dz(npoints, ratio, target_depth, min_dz=0.0001, tolerance=1):
 
 
 def angle_between(v1, v2, v3):
-    """Returns angle v2-v1-v3 that is between the vectors v1-v2 and v1-v3."""
+    """Returns the angle v2-v1-v3 (in radians). That is the angle between vectors v1-v2 and v1-v3."""
 
     v1xv2 = np.cross(v1, v2)
     v1xv3 = np.cross(v1, v3)
+
     cosangle = np.dot(v1xv2, v1xv3) / np.sqrt(
         np.dot(v1xv2, v1xv2) * np.dot(v1xv3, v1xv3)
     )

From d091f3271f88b6d3f23c5374a2d6afc4c1be8329 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 10:11:59 +1000
Subject: [PATCH 16/23] latlon_to_cartesian, quadilateral_area(s) now know R

---
 regional_mom6/regional_mom6.py | 39 ++++++++++++++++++++++------------
 1 file changed, 26 insertions(+), 13 deletions(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index 43a626d4..2edbf8ca 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -299,25 +299,36 @@ def quadilateral_area(v1, v2, v3, v4):
     orientation is implied). The area is computed as the excess of the
     sum of the spherical angles of the quadrilateral from 2π."""
 
+    if not (np.isclose(np.dot(v1, v1), np.dot(v2, v2)) &
+            np.isclose(np.dot(v1, v1), np.dot(v2, v2)) &
+            np.isclose(np.dot(v1, v1), np.dot(v3, v3)) &
+            np.isclose(np.dot(v1, v1), np.dot(v4, v4))):
+        raise Exception("vectors provided don't have same length")
+
+    R = np.sqrt(np.dot(v1, v1))
+
     a1 = angle_between(v1, v2, v4)
     a2 = angle_between(v2, v3, v1)
     a3 = angle_between(v3, v4, v2)
     a4 = angle_between(v4, v1, v3)
 
-    return a1 + a2 + a3 + a4 - 2.0 * np.pi
+    return (a1 + a2 + a3 + a4 - 2 * np.pi) * R**2
+
+
+def latlon_to_cartesian(lat, lon, R=1):
+    """Convert latitude-longitude (in degrees) to Cartesian coordinates on a sphere of radius `R`.
+    By default `R = 1`."""
 
+    x = R * np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(lon))
+    y = R * np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(lon))
+    z = R * np.sin(np.deg2rad(lat))
 
-def latlon_to_cartesian(lat, lon):
-    """Convert latitude-longitude (in degrees) to Cartesian coordinates on a unit sphere."""
-    x = np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(lon))
-    y = np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(lon))
-    z = np.sin(np.deg2rad(lat))
     return x, y, z
 
 
-def quadilateral_areas(lat, lon):
-    """Returns area of spherical quadrilaterals on the unit sphere that are formed
-    by constant latitude and longitude lines on the `lat`-`lon` grid provided.
+def quadilateral_areas(lat, lon, R=1):
+    """Returns area of spherical quadrilaterals on a sphere of radius `R`. By default, `R = 1`.
+    The quadrilaterals are formed by constant latitude and longitude lines on the `lat`-`lon` grid provided.
 
     Args:
         lat (array): Array of latitude points (in degrees)
@@ -329,7 +340,7 @@ def quadilateral_areas(lat, lon):
                        then `areas` is `(m-1) x (n-1)`.
     """
 
-    x, y, z = latlon_to_cartesian(lat, lon)
+    x, y, z = latlon_to_cartesian(lat, lon, R)
 
     nx, ny = np.shape(lat)
 
@@ -372,18 +383,20 @@ def rectangular_hgrid(λ, φ):
 
     dλ = λ[1] - λ[0]  # assuming that longitude is uniformly spaced
 
-    # dx = R * cos(φ) * np.deg2rad(dλ) / 2. Note, we divide dy by 2 because we're on the supergrid
+    # dx = R * cos(φ) * np.deg2rad(dλ) / 2
+    # Note: division by 2 because we're on the supergrid
     dx = np.broadcast_to(
         R * np.cos(np.deg2rad(φ)) * np.deg2rad(dλ) / 2,
         (λ.shape[0] - 1, φ.shape[0]),
     ).T
 
-    # dy = R * np.deg2rad(dφ) / 2. Note, we divide dy by 2 because we're on the supergrid
+    # dy = R * np.deg2rad(dφ) / 2
+    # Note: division by 2 because we're on the supergrid
     dy = np.broadcast_to(R * np.deg2rad(np.diff(φ)) / 2, (λ.shape[0], φ.shape[0] - 1)).T
 
     lon, lat = np.meshgrid(λ, φ)
 
-    area = quadilateral_areas(lat, lon) * R**2
+    area = quadilateral_areas(lat, lon, R)
 
     attrs = {
         "tile": {

From a18ff478ce16965c0c84bea1b01db18f38b09c4e Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 10:13:26 +1000
Subject: [PATCH 17/23] black formating

---
 regional_mom6/regional_mom6.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index 2edbf8ca..8529f7af 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -299,10 +299,12 @@ def quadilateral_area(v1, v2, v3, v4):
     orientation is implied). The area is computed as the excess of the
     sum of the spherical angles of the quadrilateral from 2π."""
 
-    if not (np.isclose(np.dot(v1, v1), np.dot(v2, v2)) &
-            np.isclose(np.dot(v1, v1), np.dot(v2, v2)) &
-            np.isclose(np.dot(v1, v1), np.dot(v3, v3)) &
-            np.isclose(np.dot(v1, v1), np.dot(v4, v4))):
+    if not (
+        np.isclose(np.dot(v1, v1), np.dot(v2, v2))
+        & np.isclose(np.dot(v1, v1), np.dot(v2, v2))
+        & np.isclose(np.dot(v1, v1), np.dot(v3, v3))
+        & np.isclose(np.dot(v1, v1), np.dot(v4, v4))
+    ):
         raise Exception("vectors provided don't have same length")
 
     R = np.sqrt(np.dot(v1, v1))

From 3147a518239226d7638f0966f2c28abcf1a155df Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 10:16:31 +1000
Subject: [PATCH 18/23] correct phrasing in docstring

---
 regional_mom6/regional_mom6.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index 8529f7af..91a05037 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -296,7 +296,7 @@ def angle_between(v1, v2, v3):
 def quadilateral_area(v1, v2, v3, v4):
     """Returns area of a spherical quadrilateral on the unit sphere that
     has vertices on 3-vectors `v1`, `v2`, `v3`, `v4` (counter-clockwise
-    orientation is implied). The area is computed as the excess of the
+    orientation is implied). The area is computed via the excess of the
     sum of the spherical angles of the quadrilateral from 2π."""
 
     if not (

From 139ac7092388490cb7b2a2a8922844ce364e5fe8 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 10:19:24 +1000
Subject: [PATCH 19/23] better exception error

---
 regional_mom6/regional_mom6.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index 91a05037..0533c2e9 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -305,7 +305,7 @@ def quadilateral_area(v1, v2, v3, v4):
         & np.isclose(np.dot(v1, v1), np.dot(v3, v3))
         & np.isclose(np.dot(v1, v1), np.dot(v4, v4))
     ):
-        raise Exception("vectors provided don't have same length")
+        raise Exception("vectors provided must have the same length")
 
     R = np.sqrt(np.dot(v1, v1))
 

From aa5ed7acee90ee27ef1a2f26938cff4a618fc6c0 Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 10:25:12 +1000
Subject: [PATCH 20/23] test the exception error by quadilateral_area

---
 tests/test_grid_generation.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 6dda9066..499afae0 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -57,6 +57,19 @@ def test_quadilateral_area(v1, v2, v3, v4, true_area):
     assert np.isclose(quadilateral_area(v1, v2, v3, v4), true_area)
 
 
+v1 = latlon_to_cartesian(0, 0, R=2)
+v2 = latlon_to_cartesian(90, 0, R=2)
+v3 = latlon_to_cartesian(0, 90, R=2)
+v4 = latlon_to_cartesian(-90, 0, R=2.1)
+
+
+def test_quadilateral_area_exception():
+    with pytest.raises(Exception) as excinfo:
+        quadilateral_area(v1, v2, v3, v4)
+
+    assert str(excinfo.value) == "vectors provided must have the same length"
+
+
 # create a lat-lon mesh that covers 1/4 of the North Hemisphere
 lon1, lat1 = np.meshgrid(np.linspace(0, 90, 5), np.linspace(0, 90, 5))
 area1 = 1 / 8 * (4 * np.pi)

From c423f5246c3bfff5884153b59af138c48ecb08ee Mon Sep 17 00:00:00 2001
From: "Navid C. Constantinou" <navidcy@users.noreply.github.com>
Date: Wed, 6 Sep 2023 11:12:20 +1000
Subject: [PATCH 21/23] Exception -> ValueError

---
 regional_mom6/regional_mom6.py | 2 +-
 tests/test_grid_generation.py  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index 0533c2e9..0cf55ac7 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -305,7 +305,7 @@ def quadilateral_area(v1, v2, v3, v4):
         & np.isclose(np.dot(v1, v1), np.dot(v3, v3))
         & np.isclose(np.dot(v1, v1), np.dot(v4, v4))
     ):
-        raise Exception("vectors provided must have the same length")
+        raise ValueError("vectors provided must have the same length")
 
     R = np.sqrt(np.dot(v1, v1))
 
diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 499afae0..0be9f4a8 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -64,7 +64,7 @@ def test_quadilateral_area(v1, v2, v3, v4, true_area):
 
 
 def test_quadilateral_area_exception():
-    with pytest.raises(Exception) as excinfo:
+    with pytest.raises(ValueError) as excinfo:
         quadilateral_area(v1, v2, v3, v4)
 
     assert str(excinfo.value) == "vectors provided must have the same length"

From dccdede30ddbb5396fae604ef2e9e4389f5ab52a Mon Sep 17 00:00:00 2001
From: Angus Gibson <angus@agibson.me>
Date: Wed, 6 Sep 2023 01:19:35 +0000
Subject: [PATCH 22/23] Implement quadrilateral_areas using broadcasting

We need to define a vecdot function (which is available in the Array
API but not the main numpy API, sadly), which will compute the dot
product of the last axis of two arrays.
---
 regional_mom6/regional_mom6.py | 45 +++++++++++++---------------------
 regional_mom6/utils.py         |  5 ++++
 tests/test_grid_generation.py  | 32 ++++++++++++------------
 3 files changed, 38 insertions(+), 44 deletions(-)
 create mode 100644 regional_mom6/utils.py

diff --git a/regional_mom6/regional_mom6.py b/regional_mom6/regional_mom6.py
index 0cf55ac7..a3dd6ad0 100644
--- a/regional_mom6/regional_mom6.py
+++ b/regional_mom6/regional_mom6.py
@@ -13,6 +13,7 @@
 from dask.diagnostics import ProgressBar
 import datetime as dt
 import warnings
+from .utils import vecdot
 
 warnings.filterwarnings("ignore")
 
@@ -22,8 +23,8 @@
     "dz",
     "angle_between",
     "latlon_to_cartesian",
-    "quadilateral_area",
-    "quadilateral_areas",
+    "quadrilateral_area",
+    "quadrilateral_areas",
     "rectangular_hgrid",
     "experiment",
     "segment",
@@ -286,28 +287,28 @@ def angle_between(v1, v2, v3):
     v1xv2 = np.cross(v1, v2)
     v1xv3 = np.cross(v1, v3)
 
-    cosangle = np.dot(v1xv2, v1xv3) / np.sqrt(
-        np.dot(v1xv2, v1xv2) * np.dot(v1xv3, v1xv3)
+    cosangle = vecdot(v1xv2, v1xv3) / np.sqrt(
+        vecdot(v1xv2, v1xv2) * vecdot(v1xv3, v1xv3)
     )
 
     return np.arccos(cosangle)
 
 
-def quadilateral_area(v1, v2, v3, v4):
+def quadrilateral_area(v1, v2, v3, v4):
     """Returns area of a spherical quadrilateral on the unit sphere that
     has vertices on 3-vectors `v1`, `v2`, `v3`, `v4` (counter-clockwise
     orientation is implied). The area is computed via the excess of the
     sum of the spherical angles of the quadrilateral from 2π."""
 
     if not (
-        np.isclose(np.dot(v1, v1), np.dot(v2, v2))
-        & np.isclose(np.dot(v1, v1), np.dot(v2, v2))
-        & np.isclose(np.dot(v1, v1), np.dot(v3, v3))
-        & np.isclose(np.dot(v1, v1), np.dot(v4, v4))
+        np.all(np.isclose(vecdot(v1, v1), vecdot(v2, v2)))
+        & np.all(np.isclose(vecdot(v1, v1), vecdot(v2, v2)))
+        & np.all(np.isclose(vecdot(v1, v1), vecdot(v3, v3)))
+        & np.all(np.isclose(vecdot(v1, v1), vecdot(v4, v4)))
     ):
         raise ValueError("vectors provided must have the same length")
 
-    R = np.sqrt(np.dot(v1, v1))
+    R = np.sqrt(vecdot(v1, v1))
 
     a1 = angle_between(v1, v2, v4)
     a2 = angle_between(v2, v3, v1)
@@ -328,7 +329,7 @@ def latlon_to_cartesian(lat, lon, R=1):
     return x, y, z
 
 
-def quadilateral_areas(lat, lon, R=1):
+def quadrilateral_areas(lat, lon, R=1):
     """Returns area of spherical quadrilaterals on a sphere of radius `R`. By default, `R = 1`.
     The quadrilaterals are formed by constant latitude and longitude lines on the `lat`-`lon` grid provided.
 
@@ -342,23 +343,11 @@ def quadilateral_areas(lat, lon, R=1):
                        then `areas` is `(m-1) x (n-1)`.
     """
 
-    x, y, z = latlon_to_cartesian(lat, lon, R)
+    coords = np.dstack(latlon_to_cartesian(lat, lon, R))
 
-    nx, ny = np.shape(lat)
-
-    areas = np.zeros((nx - 1, ny - 1))
-
-    for j in range(ny - 1):
-        for i in range(nx - 1):
-            # construct the 4 3-vectors v1, v2, v3, v4 that point to the vertices of the quadrilaterals
-            v1 = [x[i, j], y[i, j], z[i, j]]
-            v2 = [x[i, j + 1], y[i, j + 1], z[i, j + 1]]
-            v3 = [x[i + 1, j + 1], y[i + 1, j + 1], z[i + 1, j + 1]]
-            v4 = [x[i + 1, j], y[i + 1, j], z[i + 1, j]]
-
-            areas[i, j] = quadilateral_area(v1, v2, v3, v4)
-
-    return areas
+    return quadrilateral_area(
+        coords[:-1, :-1, :], coords[:-1, 1:, :], coords[1:, 1:, :], coords[1:, :-1, :]
+    )
 
 
 def rectangular_hgrid(λ, φ):
@@ -398,7 +387,7 @@ def rectangular_hgrid(λ, φ):
 
     lon, lat = np.meshgrid(λ, φ)
 
-    area = quadilateral_areas(lat, lon, R)
+    area = quadrilateral_areas(lat, lon, R)
 
     attrs = {
         "tile": {
diff --git a/regional_mom6/utils.py b/regional_mom6/utils.py
new file mode 100644
index 00000000..36749635
--- /dev/null
+++ b/regional_mom6/utils.py
@@ -0,0 +1,5 @@
+import numpy as np
+
+
+def vecdot(v1, v2):
+    return np.sum(v1 * v2, axis=-1)
diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index 0be9f4a8..b9148477 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -2,8 +2,8 @@
 import pytest
 from regional_mom6 import angle_between
 from regional_mom6 import latlon_to_cartesian
-from regional_mom6 import quadilateral_area
-from regional_mom6 import quadilateral_areas
+from regional_mom6 import quadrilateral_area
+from regional_mom6 import quadrilateral_areas
 from regional_mom6 import rectangular_hgrid
 import xarray as xr
 
@@ -38,23 +38,23 @@ def test_angle_between(v1, v2, v3, true_angle):
     ("v1", "v2", "v3", "v4", "true_area"),
     [
         (
-            latlon_to_cartesian(0, 0),
-            latlon_to_cartesian(0, 90),
-            latlon_to_cartesian(90, 0),
-            latlon_to_cartesian(0, -90),
+            np.dstack(latlon_to_cartesian(0, 0)),
+            np.dstack(latlon_to_cartesian(0, 90)),
+            np.dstack(latlon_to_cartesian(90, 0)),
+            np.dstack(latlon_to_cartesian(0, -90)),
             np.pi,
         ),
         (
-            latlon_to_cartesian(0, 0),
-            latlon_to_cartesian(90, 0),
-            latlon_to_cartesian(0, 90),
-            latlon_to_cartesian(-90, 0),
+            np.dstack(latlon_to_cartesian(0, 0)),
+            np.dstack(latlon_to_cartesian(90, 0)),
+            np.dstack(latlon_to_cartesian(0, 90)),
+            np.dstack(latlon_to_cartesian(-90, 0)),
             np.pi,
         ),
     ],
 )
-def test_quadilateral_area(v1, v2, v3, v4, true_area):
-    assert np.isclose(quadilateral_area(v1, v2, v3, v4), true_area)
+def test_quadrilateral_area(v1, v2, v3, v4, true_area):
+    assert np.isclose(quadrilateral_area(v1, v2, v3, v4), true_area)
 
 
 v1 = latlon_to_cartesian(0, 0, R=2)
@@ -63,9 +63,9 @@ def test_quadilateral_area(v1, v2, v3, v4, true_area):
 v4 = latlon_to_cartesian(-90, 0, R=2.1)
 
 
-def test_quadilateral_area_exception():
+def test_quadrilateral_area_exception():
     with pytest.raises(ValueError) as excinfo:
-        quadilateral_area(v1, v2, v3, v4)
+        quadrilateral_area(v1, v2, v3, v4)
 
     assert str(excinfo.value) == "vectors provided must have the same length"
 
@@ -86,8 +86,8 @@ def test_quadilateral_area_exception():
         (lat2, lon2, area2),
     ],
 )
-def test_quadilateral_areas(lat, lon, true_area):
-    assert np.isclose(np.sum(quadilateral_areas(lat, lon)), true_area)
+def test_quadrilateral_areas(lat, lon, true_area):
+    assert np.isclose(np.sum(quadrilateral_areas(lat, lon)), true_area)
 
 
 # a simple test that rectangular_hgrid runs without erroring

From 8d45e193fbc5e30f998fcf03757c6a68f6779ede Mon Sep 17 00:00:00 2001
From: Angus Gibson <angus@agibson.me>
Date: Wed, 6 Sep 2023 02:19:45 +0000
Subject: [PATCH 23/23] Fix quadrilateral area exception test

---
 tests/test_grid_generation.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/tests/test_grid_generation.py b/tests/test_grid_generation.py
index b9148477..ea0f7662 100644
--- a/tests/test_grid_generation.py
+++ b/tests/test_grid_generation.py
@@ -57,10 +57,10 @@ def test_quadrilateral_area(v1, v2, v3, v4, true_area):
     assert np.isclose(quadrilateral_area(v1, v2, v3, v4), true_area)
 
 
-v1 = latlon_to_cartesian(0, 0, R=2)
-v2 = latlon_to_cartesian(90, 0, R=2)
-v3 = latlon_to_cartesian(0, 90, R=2)
-v4 = latlon_to_cartesian(-90, 0, R=2.1)
+v1 = np.dstack(latlon_to_cartesian(0, 0, R=2))
+v2 = np.dstack(latlon_to_cartesian(90, 0, R=2))
+v3 = np.dstack(latlon_to_cartesian(0, 90, R=2))
+v4 = np.dstack(latlon_to_cartesian(-90, 0, R=2.1))
 
 
 def test_quadrilateral_area_exception():