Simplify angle_between + refactor quadrilateral_area

[navidcy]

Uses a more geometrically-inspired definition for the angle_between. Also refactors quadrilateral_area. Now quadrilateral_area computes the area of one quadrilateral and quadrilateral_areas computes the areas of all quadrilaterals on a lat-lon grid.

Closes #39

[angus-g]

Changes to angle_between look good, but not sure if you intended to use the random vector generation in the test or not?

[navidcy]

somehow the tests are failing tho... I'll figure it out

[angus-g]

Yeah, it looks like that function is actually being called with some kind of array of vectors, so it needs to be multidimension-aware?

[navidcy]

OK, I see the issue...

It's that quadilateral_area gives angle_between 3D arrays... (e.g. arrays of np.shape = (3, nx, ny)) while angle_between expects vectors of size 3.

How can I make it broadcast over the 2nd and 3rd dimension? I don't remember how to do it in python... :(

[angus-g]

I guess we could do np.cross(..., axis=0) and np.tensordot(..., axis=([0],[0])), but that's a bit ugly... It might also work if we use arrays of shape (nx, ny, 3) instead? That would also be the "correct" indexing convention for numpy -- having the vector axis first is the Fortran convention.

[navidcy]

This

    nv, nx, ny = np.shape(c0)

    c0 = np.reshape(c0, (nx, ny, nv))
    c1 = np.reshape(c1, (nx, ny, nv))
    c2 = np.reshape(c2, (nx, ny, nv))
    c3 = np.reshape(c3, (nx, ny, nv))

reshapes the inputs... But then how do I ask the angle_between to broadcast over last two dimensions?

[angus-g]

I think if you write it in such a way that the dot products and cross products are computed over the last axis, it'll just work as expected. I think that's the case for the cross product, from https://numpy.org/doc/stable/reference/generated/numpy.cross.html#numpy.cross:

The cross product of a and b in is a vector perpendicular to both a and b. If a and b are arrays of vectors, the vectors are defined by the last axis of a and b by default, and these axes can have dimensions 2 or 3.

But I'm not sure what to do for the dot product...

[navidcy]

So at the moment quadrilateral_area computes an array of areas.

I suggest we change it so that

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

that gives the area of a single quadrilateral. Then we have

def quadilateral_areas(lat, lon):
    ....

that uses quadilateral_area to computes the quadrilateral areas for the grid.

How do you like that?

[navidcy]

OK, I did it but it's not using broadcasting but it's doing the dreaded for-loop. See

regional-mom6/regional_mom6/regional_mom6.py

Lines 322 to 329 in 30254bd

	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)

@angus-g or @ashjbarnes can you convert it to broadcasting or something more efficient?

[codecov]

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.39% 🎉

Comparison is base (2f9d12b) 15.42% compared to head (8d45e19) 15.81%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #65      +/-   ##
==========================================
+ Coverage   15.42%   15.81%   +0.39%     
==========================================
  Files           2        3       +1     
  Lines         428      430       +2     
==========================================
+ Hits           66       68       +2     
  Misses        362      362              

Files Changed	Coverage Δ
regional_mom6/regional_mom6.py	14.69% <100.00%> (-0.21%)	⬇️
regional_mom6/utils.py	100.00% <100.00%> (ø)

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[angus-g]

Yeah, I'll take a look.

[navidcy]

@angus-g I added a validation in quadilateral_area via d091f32 that checks if all vectors provided are of the same length. Can we add a test to ensure that en exception is raised when vectors of different length are provided?

[navidcy]

@angus-g I added a validation in quadilateral_area via d091f32 that checks if all vectors provided are of the same length. Can we add a test to ensure that en exception is raised when vectors of different length are provided?

chatGPT told me how to do it! done via aa5ed7a

[ashjbarnes]

I like it - very verbose and readable

[navidcy]

I like it - very verbose and readable

OK, just the optimization of the for-loop with broadcasting and we are ready to merge then :)

[angus-g]

Here's an attempt at doing the broadcasting: dccdede

It's a little ugly to have to dstack everything in the tests, but on a 50x50 grid it takes 1.6ms/it vs. 715ms/it (400x speedup) vs. the for-loop.