Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 56 additions & 0 deletions src/xregrid/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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],
Expand Down
78 changes: 78 additions & 0 deletions tests/test_aero_grid_lcc.py
Original file line number Diff line number Diff line change
@@ -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
Loading