diff --git a/src/xregrid/utils.py b/src/xregrid/utils.py index 711ce2d..218cd1c 100644 --- a/src/xregrid/utils.py +++ b/src/xregrid/utils.py @@ -778,6 +778,62 @@ def create_grid_from_ioapi( return ds +def create_lcc_grid( + extent: Tuple[float, float, float, float], + res: Union[float, Tuple[float, float]], + lat_1: float, + lat_2: float, + lat_0: float, + lon_0: float, + add_bounds: bool = True, + chunks: Optional[Union[int, Dict[str, int]]] = None, +) -> xr.Dataset: + """ + Create a structured grid dataset with a Lambert Conformal Conic (LCC) projection. + + Commonly used for regional meteorological models (e.g., HRRR, NAM, CORDEX). + + Parameters + ---------- + extent : Tuple[float, float, float, float] + Grid extent in projection units (meters): (min_x, max_x, min_y, max_y). + res : float or Tuple[float, float] + Grid resolution in meters. If float, same resolution in x and y. + lat_1 : float + First standard parallel. + lat_2 : float + Second standard parallel. + lat_0 : float + Latitude of projection origin. + lon_0 : float + Longitude of projection origin (central meridian). + add_bounds : bool, default True + Whether to add cell boundary coordinates. + chunks : int or Dict[str, int], optional + Chunk sizes for the resulting dask-backed dataset. + + Returns + ------- + xr.Dataset + The grid dataset containing 'lat', 'lon' and projected coordinates 'x', 'y'. + """ + crs = ( + f"+proj=lcc +lat_1={lat_1} +lat_2={lat_2} +lat_0={lat_0} " + f"+lon_0={lon_0} +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" + ) + ds = create_grid_from_crs(crs, extent, res, add_bounds=add_bounds, chunks=chunks) + + # Backend detection for provenance + is_lazy = chunks is not None + backend = "Lazy" if is_lazy else "Eager" + + update_history( + ds, + f"Created Lambert Conformal Conic grid (lat_1={lat_1}, lat_2={lat_2}, lat_0={lat_0}, lon_0={lon_0}) using xregrid ({backend}).", + ) + return ds + + def create_sinusoidal_grid( extent: Tuple[float, float, float, float], res: Union[float, Tuple[float, float], str], diff --git a/tests/test_aero_grid_lcc.py b/tests/test_aero_grid_lcc.py new file mode 100644 index 0000000..56921b4 --- /dev/null +++ b/tests/test_aero_grid_lcc.py @@ -0,0 +1,78 @@ +import xarray as xr +import pytest +from xregrid.utils import create_lcc_grid + + +def test_lcc_grid_backend_consistency(): + """Verify that LCC grid generation yields identical results for NumPy and Dask.""" + pytest.importorskip("pyproj") + + # Define a small LCC grid + extent = (-100000, 100000, -100000, 100000) + res = 20000 + lat_1, lat_2 = 33, 45 + lat_0, lon_0 = 40, -97 + + ds_eager = create_lcc_grid( + extent=extent, + res=res, + lat_1=lat_1, + lat_2=lat_2, + lat_0=lat_0, + lon_0=lon_0, + chunks=None, + ) + + ds_lazy = create_lcc_grid( + extent=extent, + res=res, + lat_1=lat_1, + lat_2=lat_2, + lat_0=lat_0, + lon_0=lon_0, + chunks={"x": 5, "y": 5}, + ) + + # 1. Numerical Consistency + xr.testing.assert_allclose(ds_eager, ds_lazy.compute()) + + # 2. Laziness Check + # In xarray, dimension coordinates are eager. + # Non-dimension coordinates and DataArrays should be lazy. + assert hasattr(ds_lazy.lat.data, "dask") + assert hasattr(ds_lazy.lon.data, "dask") + assert hasattr(ds_lazy.lat_b.data, "dask") + assert hasattr(ds_lazy.lon_b.data, "dask") + + # 3. Metadata and CF-compliance + assert "lat" in ds_eager.coords + assert "lon" in ds_eager.coords + assert ds_eager.lat.attrs["standard_name"] == "latitude" + assert ds_eager.lon.attrs["standard_name"] == "longitude" + assert "crs" in ds_eager.attrs + assert "Lambert Conic Conformal" in ds_eager.attrs["crs"] + + # 4. Provenance + assert "history" in ds_eager.attrs + assert "Created Lambert Conformal Conic grid" in ds_eager.attrs["history"] + assert "(Eager)" in ds_eager.attrs["history"] + assert "(Lazy)" in ds_lazy.attrs["history"] + + +def test_lcc_grid_resolution_tuple(): + """Verify LCC grid works with a tuple for resolution.""" + pytest.importorskip("pyproj") + extent = (-10000, 10000, -10000, 10000) + res = (1000, 2000) + + ds = create_lcc_grid( + extent=extent, + res=res, + lat_1=30, + lat_2=60, + lat_0=45, + lon_0=-100, + ) + + assert ds.x.size == 20 # 20000 / 1000 + assert ds.y.size == 10 # 20000 / 2000