angle_between and quad_area functions

[navidcy]

I'm wondering if we really need to write our own functions for angle between vectors and the surface of a spherical quadrilateral. I have the gut feeling that there should be a python package that does that. Then we don't have to maintain and test these functions:

https://github.com/COSIMA/mom6-regional/blob/145e9f2b75b1009c70c4d00a4b32d14714870bf6/mom6_regional/mom6_regional.py#L264-L298

If we really need to rewrite those then a docstring is needed. ;)

[navidcy]

Regarding the algorithm that computes the angle between vectors, I'd suggest

v1xv2 = np.cross(v1, v2)
v1xv3 = np.cross(v1, v3)
cosangle = np.dot(v1xv2, v1xv3) / np.sqrt( np.dot(v1xv2, v1xv1) * np.dot(v1xv3, v1xv3) )
angle = np.arccos(cosangle)

This way we reuse numpy's cross and dot functions.

[ashjbarnes]

Yeah good points Navid! That would be ideal. I just copied and pasted for now to be speedy and not get bogged down in details but this is certainly more elegant.

[navidcy]

Fair enough!

OK. Well we can first add tests for function as they are currently implemented and then changing will be trivial (if we break anything tests will fail).

[navidcy]

Btw, as is defined at the moment, angle_between fails if any of the vectors is the origin [0, 0, 0] (everything is zero and we end up dividing with 0).

[ashjbarnes]

Golden rule of (Ashley’s) regional modelling: no poles

…

On Thu, 22 Jun 2023 at 18:11, Navid C. Constantinou < ***@***.***> wrote: Btw, as is defined at the moment, angle_between fails if any of the vectors is the origin [0, 0, 0] (everything is zero and we end up dividing with 0). — Reply to this email directly, view it on GitHub <#39 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AMWQL4C3E3BK5RDQU3S4UFTXMP4Z7ANCNFSM6AAAAAAZHEVUBE> . You are receiving this because you commented.Message ID: ***@***.***>

[navidcy]

hm... (0, 0, 0) is the Earth's center, right?

[ashjbarnes]

Oh true! Yeah definitely a no go zone I think

…

On Thu, 22 Jun 2023 at 18:55, Navid C. Constantinou < ***@***.***> wrote: hm... (0, 0, 0) is the Earth's center, right? — Reply to this email directly, view it on GitHub <#39 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AMWQL4ACHHFFIIL3FTD7SVLXMQB7HANCNFSM6AAAAAAZHEVUBE> . You are receiving this because you commented.Message ID: ***@***.***>

[navidcy]

It seems that the angle_between is only used by quad_area, right?

And that the quad_area is a function that computes the horizontal area of a cell centered at (lat, lon). The shape of that cell is a spherical quadrilateral.

Perhaps we can utilize exact formulas for the area of spherical triangles. I recently coded these for the spherical grid in Oceananigans:
https://github.com/CliMA/Oceananigans.jl/blob/eb38eeade577eca5056b99ca8839ca6c674ae0e2/src/Grids/grid_utils.jl#L548-L664

The relevant method of the above for the area of a spherical triangle is the one that returns the area of the triangle given the three 3-vectors of its vertices:
https://github.com/CliMA/Oceananigans.jl/blob/eb38eeade577eca5056b99ca8839ca6c674ae0e2/src/Grids/grid_utils.jl#L594

The area of the quadrilateral is then computed as 1/2 x sum of area of all possible combinations of triangles you can form with the 4 vertices. Normally, the quadrilateral is just the sum of two of the triangles. But by summing all 4 possible triangles and then taking 1/2 of that area we bypass the need to ensure that we selected the vertices to form the triangles in a way that the triangles don't overlap.

[navidcy]

If I understand correct what's happening here, quad_area computes the area of the spherical quadrilateral based on the formula

area = A + B + C + D - 2π

where A,... are the angles on the vertices of the quadrilateral.

BUT the angles then are computed via an approximation? Connect the consecutive points with an arrow and compute the angle between those arrows? This is only approximately correct as the resolution increases, but for coarse resolution there is some error, right?

A good test is to define a grid that covers a know portion of the globe and ensure that the sum of all the areas is exactly that. E.g., λ ∈ [-45ᵒ, 45ᵒ] and φ ∈ [0, 30ᵒ] should give an area of πR²/4.

[ashjbarnes]

Was this resolved by this commit?

[navidcy]

Nope. I still this can be simplified as in #39 (comment)

Now with the tests perhaps we can give it a go.