|
| 1 | +""" |
| 2 | +mom6_static.py: Utilities for reading and extracting grid information from MOM6 static files. |
| 3 | +""" |
| 4 | + |
| 5 | +import xarray as xr |
| 6 | +import numpy as np |
| 7 | + |
| 8 | + |
| 9 | +def compute_cell_bounds_from_corners(corner_array): |
| 10 | + """ |
| 11 | + Given a 2D array of cell corners (shape [ny+1, nx+1]), |
| 12 | + compute bounds for each cell as [min, max]. |
| 13 | + Returns a (ny*nx, 2) array for CMOR. |
| 14 | + """ |
| 15 | + # For each cell, get the 4 corners and compute min/max |
| 16 | + ny, nx = corner_array.shape[0] - 1, corner_array.shape[1] - 1 |
| 17 | + print(f"nx = {nx}, ny = {ny}") |
| 18 | + bounds = np.empty((ny * nx, 2), dtype=corner_array.dtype) |
| 19 | + idx = 0 |
| 20 | + for j in range(ny): |
| 21 | + for i in range(nx): |
| 22 | + cell_corners = [ |
| 23 | + corner_array[j, i], |
| 24 | + corner_array[j, i + 1], |
| 25 | + corner_array[j + 1, i], |
| 26 | + corner_array[j + 1, i + 1], |
| 27 | + ] |
| 28 | + bounds[idx, 0] = np.nanmin(cell_corners) |
| 29 | + bounds[idx, 1] = np.nanmax(cell_corners) |
| 30 | + idx += 1 |
| 31 | + print(f"[MOM6 bounds debug] ny: {ny}, nx: {nx}, bounds shape: {bounds.shape}") |
| 32 | + return bounds |
| 33 | + |
| 34 | + |
| 35 | +def load_mom6_static(static_path): |
| 36 | + """Load MOM6 static file and return geolat, geolon, geolat_c, geolon_c arrays.""" |
| 37 | + ds = xr.open_dataset(static_path) |
| 38 | + # For a supergrid, centers are every other point, bounds are full array |
| 39 | + # Target: centers (480, 540), corners (481, 541) |
| 40 | + # For a supergrid, centers are every other point starting at 1, |
| 41 | + # bounds are every other point starting at 0 |
| 42 | + geolat = ds["y"].values[1::2, 1::2] # (480, 540) |
| 43 | + geolon = ds["x"].values[1::2, 1::2] # (480, 540) |
| 44 | + geolat_c = ds["y"].values[0::2, 0::2] # (481, 541) |
| 45 | + geolon_c = ds["x"].values[0::2, 0::2] # (481, 541) |
| 46 | + # Normalize longitudes to 0-360 |
| 47 | + geolon = np.mod(geolon, 360.0) |
| 48 | + geolon_c = np.mod(geolon_c, 360.0) |
| 49 | + # Compute bounds arrays using corners |
| 50 | + lat_bnds = compute_cell_bounds_from_corners(geolat_c) |
| 51 | + lon_bnds = compute_cell_bounds_from_corners(geolon_c) |
| 52 | + # Debug: print min, max, and check monotonicity |
| 53 | + print( |
| 54 | + f"lat_bnds shape: {lat_bnds.shape}, min: {np.nanmin(lat_bnds)}, max: {np.nanmax(lat_bnds)}" |
| 55 | + ) |
| 56 | + print( |
| 57 | + f"lon_bnds shape: {lon_bnds.shape}, min: {np.nanmin(lon_bnds)}, max: {np.nanmax(lon_bnds)}" |
| 58 | + ) |
| 59 | + return geolat, geolon, geolat_c, geolon_c |
0 commit comments