Start of Angle Calc

regional_mom6/regridding.py

@@ -379,3 +379,82 @@ def generate_encoding(
             }
 
     return encoding_dict
+
+
+
+def modulo_around_point(x, xc, Lx):
+    """
+    This function calculates the modulo around a point. Return the modulo value of x in an interval [xc-(Lx/2) xc+(Lx/2)]. If Lx<=0, then it returns x without applying modulo arithmetic.
+    Parameters
+    ----------
+    x: float
+        Value to which to apply modulo arithmetic
+    xc: float
+        Center of modulo range
+    Lx: float
+        Modulo range width
+    Returns
+    -------
+    float
+        x shifted by an integer multiple of Lx to be close to xc, 
+    """
+    if Lx <= 0:
+        return x
+    else:
+        return ((x - (xc - 0.5*Lx)) % Lx )- Lx/2 + xc
+
+
+def initialize_grid_rotation_angle(hgrid: xr.Dataset) -> xr.Dataset:
+    """
+    Calculate the angle_dx in degrees from the true x (east?) direction counterclockwise) and save as angle_dx_mom6
+    Parameters
+    ----------
+    hgrid: xr.Dataset
+        The hgrid dataset
+    Returns
+    -------
+    xr.Dataset
+        The dataset with the mom6_angle_dx variable added
+    """
+    regridding_logger.info("Initializing grid rotation angle")
+    # Direct Translation
+    pi_720deg = np.arctan(1)/180 # One quarter the conversion factor from degrees to radians
+
+    ## Check length of longitude
+    len_lon = 360.0
+    G_len_lon = hgrid.x.max() - hgrid.x.min()
+    if G_len_lon != 360:
+        regridding_logger.info("This is a regional case")
+        len_lon = G_len_lon
+
+    angles_arr = np.zeros((len(hgrid.nyp), len(hgrid.nxp)))
+
+    # Compute lonB for all points
+    lonB = np.zeros((2, 2, len(hgrid.nyp)-2, len(hgrid.nxp)-2))
+
+    # Vectorized Compute lonB - Fortran is 1-indexed, so we need to subtract 1 from the indices in lonB[m-1,n-1]
+    for n in np.arange(1,3):
+        for m in np.arange(1,3):
+            lonB[m-1, n-1] = modulo_around_point(hgrid.x[np.arange((m-2+1),(m-2+len(hgrid.nyp)-1)), np.arange((n-2+1),(n-2+len(hgrid.nxp)-1))], hgrid.x[1:-1,1:-1], len_lon)
+
+    # Vectorized Compute lon_scale
+    lon_scale = np.cos(pi_720deg* ((hgrid.y[0:-2, 0:-2] + hgrid.y[1:-1, 1:-1]) + (hgrid.y[1:-1, 0:-2] + hgrid.y[0:-2, 1:-1])))
+
+    # Vectorized Compute angle
+    angle = np.arctan2(
+        lon_scale * ((lonB[0, 1] - lonB[1, 0]) + (lonB[1, 1] - lonB[0, 0])),
+        (hgrid.y[0:-2, 1:-1] - hgrid.y[1:-1, 0:-2]) + (hgrid.y[1:-1, 1:-1] - hgrid.y[0:-2, 0:-2])
+    )
+
+    # Assign angle to angles_arr
+    angles_arr[1:-1,1:-1] = 90 - np.rad2deg(angle)
+
+
+    # Assign angles_arr to hgrid
+    hgrid["angle_dx_mom6"] = (("nyp", "nxp"), angles_arr)
+    hgrid["angle_dx_mom6"].attrs["_FillValue"] = np.nan
+    hgrid["angle_dx_mom6"].attrs["units"] = "deg"
+    hgrid["angle_dx_mom6"].attrs["description"] = "MOM6 calculates angles internally, this field replicates that for rotating boundary conditions. Use this over other angle fields for MOM6 applications"
+
+    return hgrid
+

regional_mom6/testing_to_be_deleted/angle_calc.md

@@ -0,0 +1,11 @@
+# MOM6 Angle Calculation Steps
+
+1. Calculate pi/4rads / 180 degress  = Gives a 1/4 conversion of degrees to radians. I.E. multiplying an angle in degrees by this gives the conversion to radians at 1/4 the value. 
+2. Figure out the longitudunal extent of our domain, or periodic range of longitudes. For global cases it is len_lon = 360, for our regional cases it is given by the hgrid.
+3. At each point on our hgrid, we find the point to the left, bottom left diag, bottom, and itself. We adjust each of these longitudes to be in the range of len_lon around the point itself. (module_around_point)
+4. We then find the lon_scale, which is the "trigonometric scaling factor converting changes in longitude to equivalent distances in latitudes". Whatever that actually means is we add the latitude of all four of these points from part 3 and basically average it and convert to radians. We then take the cosine of it. 
+5. Then we calculate the angle. This is a simple arctan2 so y/x. 
+    1. The "y" component is the addition of the difference between the diagonals in longitude of lonB multiplied by the lon_scale, which is our conversion to latitude.
+    2. The "x" component is the same addition of differences in latitude.
+    3. Thus, given the same units, we can call arctan to get the angle in degrees
+6. Challenge: Because this takes the left & bottom points, we can't calculate the angle at the left and bottom edges. Therefore, we can always calculate it the other way by using the right and top points.