Utilities¶

def create_global_grid(
    res_lat: float,
    res_lon: float,
    add_bounds: bool = True,
    chunks: Optional[Union[int, Dict[str, int]]] = None,
) -> xr.Dataset:
    """
    Create a global rectilinear grid dataset.

    Parameters
    ----------
    res_lat : float
        Latitude resolution in degrees.
    res_lon : float
        Longitude resolution in degrees.
    add_bounds : bool, default True
        Whether to add cell boundary coordinates.
    chunks : int or dict, optional
        Chunk sizes for the resulting dask-backed dataset.
        If None (default), returns an eager NumPy-backed dataset.

    Returns
    -------
    xr.Dataset
        The global grid dataset containing 'lat' and 'lon'.
    """
    return _create_rectilinear_grid(
        lat_range=(-90, 90),
        lon_range=(0, 360),
        res_lat=res_lat,
        res_lon=res_lon,
        add_bounds=add_bounds,
        chunks=chunks,
        history_msg=f"Created global grid ({res_lat}x{res_lon}) using xregrid.",
    )

def create_regional_grid(
    lat_range: Tuple[float, float],
    lon_range: Tuple[float, float],
    res_lat: float,
    res_lon: float,
    add_bounds: bool = True,
    chunks: Optional[Union[int, Dict[str, int]]] = None,
) -> xr.Dataset:
    """
    Create a regional rectilinear grid dataset.

    Parameters
    ----------
    lat_range : tuple of float
        (min_lat, max_lat).
    lon_range : tuple of float
        (min_lon, max_lon).
    res_lat : float
        Latitude resolution in degrees.
    res_lon : float
        Longitude resolution in degrees.
    add_bounds : bool, default True
        Whether to add cell boundary coordinates.
    chunks : int or dict, optional
        Chunk sizes for the resulting dask-backed dataset.
        If None (default), returns an eager NumPy-backed dataset.

    Returns
    -------
    xr.Dataset
        The regional grid dataset containing 'lat' and 'lon'.
    """
    return _create_rectilinear_grid(
        lat_range=lat_range,
        lon_range=lon_range,
        res_lat=res_lat,
        res_lon=res_lon,
        add_bounds=add_bounds,
        chunks=chunks,
        history_msg=f"Created regional grid ({res_lat}x{res_lon}) using xregrid.",
    )

def create_grid_like(
    obj: Union[xr.DataArray, xr.Dataset],
    res: Union[float, Tuple[float, float]],
    add_bounds: bool = True,
    chunks: Optional[Union[int, Dict[str, int]]] = None,
    extent: Optional[Tuple[float, float, float, float]] = None,
    crs: Optional[Union[str, int, Any]] = None,
) -> xr.Dataset:
    """
    Create a new grid dataset with the same extent and CRS as an existing object.

    Automatically detects the CRS and spatial extent of the input object.
    Supports both geographic (lat-lon) and projected coordinate systems.
    Efficiently handles Dask-backed objects by batching metadata discovery.

    Parameters
    ----------
    obj : xr.DataArray or xr.Dataset
        The input object to use as a template.
    res : Union[float, Tuple[float, float]]
        New grid resolution in the coordinate system units.
        If tuple, (res_x, res_y) or (res_lon, res_lat).
    add_bounds : bool, default True
        Whether to add cell boundary coordinates.
    chunks : Union[int, Dict[str, int]], optional
        Chunk sizes for the resulting dask-backed dataset.
    extent : Tuple[float, float, float, float], optional
        Override the detected extent (min_x, max_x, min_y, max_y).
        Use this to avoid hidden dask.compute() if you already know the extent.
    crs : Union[str, int, pyproj.CRS], optional
        Override the detected CRS.

    Returns
    -------
    xr.Dataset
        The new grid dataset with consistent metadata.
    """
    if crs is not None:
        if pyproj is not None:
            crs_obj = pyproj.CRS(crs)
        else:
            crs_obj = crs
    else:
        crs_obj = get_crs_info(obj)

    if isinstance(res, (int, float)):
        res_x = res_y = float(res)
    else:
        res_x, res_y = map(float, res)

    obj_name = getattr(obj, "name", "input")
    history_msg_base = f"Created grid like {obj_name} using xregrid."
    if hasattr(obj, "attrs") and "history" in obj.attrs:
        history_msg_base += f"\nTemplate history:\n{obj.attrs['history']}"

    if extent is not None:
        if crs_obj is None or (
            hasattr(crs_obj, "is_geographic") and crs_obj.is_geographic
        ):
            # Lat-Lon
            return _create_rectilinear_grid(
                (extent[2], extent[3]),  # lat_range
                (extent[0], extent[1]),  # lon_range
                res_y,  # res_lat
                res_x,  # res_lon
                add_bounds=add_bounds,
                chunks=chunks,
                crs=crs_obj.to_wkt() if hasattr(crs_obj, "to_wkt") else "EPSG:4326",
                history_msg=history_msg_base + " (Override Extent).",
            )
        else:
            # Projected
            return create_grid_from_crs(
                crs_obj, extent, (res_x, res_y), add_bounds=add_bounds, chunks=chunks
            )

    # Aero Optimization: Try to find extent in metadata before falling back to compute.
    if extent is None:
        try:
            # Check for geospatial metadata (e.g. from NetCDF or GeoTIFF)
            if "geospatial_lat_min" in obj.attrs and "geospatial_lat_max" in obj.attrs:
                lat_range = (
                    float(obj.attrs["geospatial_lat_min"]),
                    float(obj.attrs["geospatial_lat_max"]),
                )
                lon_range = (
                    float(obj.attrs.get("geospatial_lon_min", 0)),
                    float(obj.attrs.get("geospatial_lon_max", 360)),
                )
                extent = (lon_range[0], lon_range[1], lat_range[0], lat_range[1])
        except (KeyError, ValueError):
            pass

    if extent is not None:
        if crs_obj is None or (
            hasattr(crs_obj, "is_geographic") and crs_obj.is_geographic
        ):
            # Lat-Lon
            return _create_rectilinear_grid(
                (extent[2], extent[3]),  # lat_range
                (extent[0], extent[1]),  # lon_range
                res_y,  # res_lat
                res_x,  # res_lon
                add_bounds=add_bounds,
                chunks=chunks,
                crs=crs_obj.to_wkt() if hasattr(crs_obj, "to_wkt") else "EPSG:4326",
                history_msg=history_msg_base + " (Override Extent).",
            )
        else:
            # Projected
            return create_grid_from_crs(
                crs_obj, extent, (res_x, res_y), add_bounds=add_bounds, chunks=chunks
            )

    # Discovery logic: we need min/max. We use batch compute if lazy to minimize roundtrips.
    # Aero-Optimization: Use Xarray indexes for 1D dimension coordinates to avoid hidden computes.

    # 1. Try to find projected coordinates
    try:
        x_da = obj.cf["projection_x_coordinate"]
        y_da = obj.cf["projection_y_coordinate"]

        try:
            # Use bounds for exact extent if available
            x_b = obj.cf.get_bounds("projection_x_coordinate")
            y_b = obj.cf.get_bounds("projection_y_coordinate")

            x_min, x_max, x_eager = _get_min_max_lazy_aware(x_b)
            y_min, y_max, y_eager = _get_min_max_lazy_aware(y_b)

            if not (x_eager and y_eager):
                if is_lazy(obj):
                    warnings.warn(
                        f"Triggering hidden compute in create_grid_like to determine extent of {obj_name}. "
                        "To avoid this, provide 'extent' explicitly."
                    )
                tasks_dict = {}
                if not x_eager:
                    tasks_dict["xmin"] = x_min
                    tasks_dict["xmax"] = x_max
                if not y_eager:
                    tasks_dict["ymin"] = y_min
                    tasks_dict["ymax"] = y_max

                results = dask.compute(tasks_dict)[0]
                extent = (
                    float(results.get("xmin", x_min)),
                    float(results.get("xmax", x_max)),
                    float(results.get("ymin", y_min)),
                    float(results.get("ymax", y_max)),
                )
            else:
                extent = (float(x_min), float(x_max), float(y_min), float(y_max))
        except Exception:
            # Fallback to centers
            x_min, x_max, x_eager = _get_min_max_lazy_aware(x_da)
            y_min, y_max, y_eager = _get_min_max_lazy_aware(y_da)

            if not (x_eager and y_eager):
                if is_lazy(obj):
                    warnings.warn(
                        f"Triggering hidden compute in create_grid_like to determine extent of {obj_name}. "
                        "To avoid this, provide 'extent' explicitly."
                    )
                # Batch everything!
                tasks_dict = {}
                if not x_eager:
                    tasks_dict["x_min"] = x_min
                    tasks_dict["x_max"] = x_max
                if not y_eager:
                    tasks_dict["y_min"] = y_min
                    tasks_dict["y_max"] = y_max

                if x_da.size > 1:
                    tasks_dict["res_x"] = abs(x_da.diff(x_da.dims[0]).mean())
                if y_da.size > 1:
                    tasks_dict["res_y"] = abs(y_da.diff(y_da.dims[0]).mean())

                results = dask.compute(tasks_dict)[0]
                x_min_val = float(results.get("x_min", x_min))
                x_max_val = float(results.get("x_max", x_max))
                y_min_val = float(results.get("y_min", y_min))
                y_max_val = float(results.get("y_max", y_max))

                res_x_orig = float(results.get("res_x", 0))
                res_y_orig = float(
                    results.get("res_y", res_x_orig if res_x_orig else 0)
                )

                extent = (
                    x_min_val - res_x_orig / 2,
                    x_max_val + res_x_orig / 2,
                    y_min_val - res_y_orig / 2,
                    y_max_val + res_y_orig / 2,
                )
            else:
                x_min_val, x_max_val = float(x_min), float(x_max)
                y_min_val, y_max_val = float(y_min), float(y_max)

                res_x_orig = (
                    abs(float(x_da.diff(x_da.dims[0]).mean())) if x_da.size > 1 else 0
                )
                res_y_orig = (
                    abs(float(y_da.diff(y_da.dims[0]).mean()))
                    if y_da.size > 1
                    else res_x_orig
                )
                extent = (
                    x_min_val - res_x_orig / 2,
                    x_max_val + res_x_orig / 2,
                    y_min_val - res_y_orig / 2,
                    y_max_val + res_y_orig / 2,
                )

        if crs_obj is None:
            # Fallback to generic geographic if no CRS found
            crs_obj = "EPSG:4326"

        return create_grid_from_crs(
            crs_obj, extent, (res_x, res_y), add_bounds=add_bounds, chunks=chunks
        )

    except (KeyError, AttributeError, ValueError):
        pass

    # 2. Fallback to Geographic (Lat-Lon)
    try:
        lat_da = _find_coord(obj, "latitude")
        lon_da = _find_coord(obj, "longitude")
        if lat_da is None or lon_da is None:
            raise KeyError("Coordinates not found")

        try:
            lat_b = obj.cf.get_bounds("latitude")
            lon_b = obj.cf.get_bounds("longitude")

            lat_min, lat_max, lat_eager = _get_min_max_lazy_aware(lat_b)
            lon_min, lon_max, lon_eager = _get_min_max_lazy_aware(lon_b)

            if not (lat_eager and lon_eager):
                if is_lazy(obj):
                    warnings.warn(
                        f"Triggering hidden compute in create_grid_like to determine extent of {obj_name}. "
                        "To avoid this, provide 'extent' explicitly."
                    )
                tasks_dict = {}
                if not lat_eager:
                    tasks_dict["lat_min"] = lat_min
                    tasks_dict["lat_max"] = lat_max
                if not lon_eager:
                    tasks_dict["lon_min"] = lon_min
                    tasks_dict["lon_max"] = lon_max

                results = dask.compute(tasks_dict)[0]
                lat_range = (
                    float(results.get("lat_min", lat_min)),
                    float(results.get("lat_max", lat_max)),
                )
                lon_range = (
                    float(results.get("lon_min", lon_min)),
                    float(results.get("lon_max", lon_max)),
                )
            else:
                lat_range = (float(lat_min), float(lat_max))
                lon_range = (float(lon_min), float(lon_max))
        except Exception:
            # Heuristic for resolution to calculate extent from centers
            lat_min, lat_max, lat_eager = _get_min_max_lazy_aware(lat_da)
            lon_min, lon_max, lon_eager = _get_min_max_lazy_aware(lon_da)

            if not (lat_eager and lon_eager):
                if is_lazy(obj):
                    warnings.warn(
                        f"Triggering hidden compute in create_grid_like to determine extent of {obj_name}. "
                        "To avoid this, provide 'extent' explicitly."
                    )
                tasks_dict = {}
                if not lat_eager:
                    tasks_dict["lat_min"] = lat_min
                    tasks_dict["lat_max"] = lat_max
                if not lon_eager:
                    tasks_dict["lon_min"] = lon_min
                    tasks_dict["lon_max"] = lon_max

                if lat_da.size > 1:
                    tasks_dict["res_lat"] = abs(lat_da.diff(lat_da.dims[0]).mean())
                if lon_da.size > 1:
                    tasks_dict["res_lon"] = abs(lon_da.diff(lon_da.dims[-1]).mean())

                results = dask.compute(tasks_dict)[0]
                lat_min_val = float(results.get("lat_min", lat_min))
                lat_max_val = float(results.get("lat_max", lat_max))
                lon_min_val = float(results.get("lon_min", lon_min))
                lon_max_val = float(results.get("lon_max", lon_max))

                res_lat_orig = float(results.get("res_lat", 0))
                res_lon_orig = float(
                    results.get("res_lon", res_lat_orig if res_lat_orig else 0)
                )

                lat_range = (
                    lat_min_val - res_lat_orig / 2,
                    lat_max_val + res_lat_orig / 2,
                )
                lon_range = (
                    lon_min_val - res_lon_orig / 2,
                    lon_max_val + res_lon_orig / 2,
                )
            else:
                lat_min_val, lat_max_val = float(lat_min), float(lat_max)
                lon_min_val, lon_max_val = float(lon_min), float(lon_max)

                res_lat_orig = (
                    abs(float(lat_da.diff(lat_da.dims[0]).mean()))
                    if lat_da.size > 1
                    else 0
                )
                res_lon_orig = (
                    abs(float(lon_da.diff(lon_da.dims[-1]).mean()))
                    if lon_da.size > 1
                    else res_lat_orig
                )
                lat_range = (
                    lat_min_val - res_lat_orig / 2,
                    lat_max_val + res_lat_orig / 2,
                )
                lon_range = (
                    lon_min_val - res_lon_orig / 2,
                    lon_max_val + res_lon_orig / 2,
                )

        return _create_rectilinear_grid(
            lat_range,
            lon_range,
            res_y,  # res_lat
            res_x,  # res_lon
            add_bounds=add_bounds,
            chunks=chunks,
            crs=crs_obj.to_wkt() if crs_obj else "EPSG:4326",
            history_msg=history_msg_base,
        )
    except (KeyError, AttributeError, ValueError):
        raise ValueError(
            "Could not detect spatial coordinates (latitude/longitude or "
            "projection_x/y) in input object."
        )

def create_grid_from_crs(
    crs: Union[str, int, Any],
    extent: Tuple[float, float, float, float],
    res: Union[float, Tuple[float, float]],
    add_bounds: bool = True,
    chunks: Optional[Union[int, Dict[str, int]]] = None,
) -> xr.Dataset:
    """
    Create a structured grid dataset from a CRS and extent.

    Parameters
    ----------
    crs : str, int, or pyproj.CRS
        The CRS of the grid (Proj4 string, EPSG code, WKT, or CRS object).
    extent : Tuple[float, float, float, float]
        Grid extent in CRS units: (min_x, max_x, min_y, max_y).
    res : float or Tuple[float, float]
        Grid resolution in CRS units. If float, same resolution in x and y.
        If tuple, (res_x, res_y).
    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'.
    """
    if isinstance(res, (int, float)):
        res_x = res_y = float(res)
    else:
        res_x, res_y = map(float, res)

    x_chunks = chunks.get("x", -1) if isinstance(chunks, dict) else chunks
    y_chunks = chunks.get("y", -1) if isinstance(chunks, dict) else chunks

    # Generate 1D coordinates in projected space
    x = _lazy_arange(extent[0] + res_x / 2, extent[1], res_x, chunks=x_chunks)
    y = _lazy_arange(extent[2] + res_y / 2, extent[3], res_y, chunks=y_chunks)

    x_da = xr.DataArray(x, dims=["x"], name="x")
    y_da = xr.DataArray(y, dims=["y"], name="y")

    # Use xr.broadcast for lazy 2D arrays
    yy_da, xx_da = xr.broadcast(y_da, x_da)

    # Ensure (y, x) order
    yy_da = yy_da.transpose("y", "x")
    xx_da = xx_da.transpose("y", "x")

    if pyproj is None:
        raise ImportError("pyproj is required for create_grid_from_crs.")
    crs_obj = pyproj.CRS(crs)

    lon, lat = xr.apply_ufunc(
        _transform_coords,
        xx_da,
        yy_da,
        kwargs={"crs_in": crs_obj},
        dask="parallelized",
        output_dtypes=[float, float],
        input_core_dims=[[], []],
        output_core_dims=[[], []],
    )

    try:
        units = crs_obj.axis_info[0].unit_name or "m"
    except (IndexError, AttributeError):
        units = "m"

    ds = xr.Dataset(
        coords={
            "y": (
                ["y"],
                y,
                {"units": units, "standard_name": "projection_y_coordinate"},
            ),
            "x": (
                ["x"],
                x,
                {"units": units, "standard_name": "projection_x_coordinate"},
            ),
            "lat": (
                ["y", "x"],
                lat.data,
                {"units": "degrees_north", "standard_name": "latitude"},
            ),
            "lon": (
                ["y", "x"],
                lon.data,
                {"units": "degrees_east", "standard_name": "longitude"},
            ),
        }
    )

    ds.attrs["crs"] = crs_obj.to_wkt()

    if add_bounds:
        xp = _get_array_namespace(x, y)
        x_b_raw = xp.stack([x - res_x / 2, x + res_x / 2, x + res_x / 2, x - res_x / 2])
        y_b_raw = xp.stack([y - res_y / 2, y - res_y / 2, y + res_y / 2, y + res_y / 2])

        x_b_da = xr.DataArray(x_b_raw, dims=["nv", "x"])
        y_b_da = xr.DataArray(y_b_raw, dims=["nv", "y"])

        yy_b_da, xx_b_da = xr.broadcast(y_b_da, x_b_da)

        lon_b, lat_b = xr.apply_ufunc(
            _transform_coords,
            xx_b_da,
            yy_b_da,
            kwargs={"crs_in": crs_obj},
            dask="parallelized",
            output_dtypes=[float, float],
            input_core_dims=[[], []],
            output_core_dims=[[], []],
        )

        ds.coords["lat_b"] = (
            ["y", "x", "nv"],
            lat_b.data.transpose(1, 2, 0),
            {"units": "degrees_north"},
        )
        ds.coords["lon_b"] = (
            ["y", "x", "nv"],
            lon_b.data.transpose(1, 2, 0),
            {"units": "degrees_east"},
        )
        ds["lat"].attrs["bounds"] = "lat_b"
        ds["lon"].attrs["bounds"] = "lon_b"

        # Add 1D projected bounds using backend-agnostic xarray operations
        x_da_1d = xr.DataArray(x, dims=["x"])
        y_da_1d = xr.DataArray(y, dims=["y"])

        # Create (N, 2) bounds
        x_b_1d = xr.concat(
            [x_da_1d - res_x / 2, x_da_1d + res_x / 2], dim="nbounds"
        ).transpose("x", "nbounds")
        y_b_1d = xr.concat(
            [y_da_1d - res_y / 2, y_da_1d + res_y / 2], dim="nbounds"
        ).transpose("y", "nbounds")

        ds.coords["x_b"] = (["x", "nbounds"], x_b_1d.data, {"units": units})
        ds.coords["y_b"] = (["y", "nbounds"], y_b_1d.data, {"units": units})
        ds["x"].attrs["bounds"] = "x_b"
        ds["y"].attrs["bounds"] = "y_b"

    # Backend detection for provenance
    is_lazy = chunks is not None
    backend = "Lazy" if is_lazy else "Eager"

    # Add extra metadata about the generated bounds if present
    bounds_msg = " with cell bounds" if add_bounds else ""
    update_history(
        ds, f"Created grid from CRS {crs} using xregrid ({backend}){bounds_msg}."
    )
    if chunks is not None:
        ds = ds.chunk(chunks)
    return ds

def create_grid_from_ioapi(
    metadata: Dict[str, Any],
    add_bounds: bool = True,
    chunks: Optional[Union[int, Dict[str, int]]] = None,
) -> xr.Dataset:
    """
    Create a structured grid dataset from IOAPI-compliant metadata.

    Supports GDTYP:
    - 1: Lat-Lon
    - 2: Lambert Conformal
    - 3: Mercator
    - 4: Stereographic
    - 5: UTM
    - 6: Polar Stereographic
    - 7: Equatorial Mercator
    - 8: Transverse Mercator
    - 9: Albers Equal Area
    - 10: Lambert Azimuthal Equal Area
    - 13: Sinusoidal

    Parameters
    ----------
    metadata : dict
        IOAPI metadata containing GDTYP, P_ALP, P_BET, P_GAM, XCENT, YCENT,
        XORIG, YORIG, XCELL, YCELL, NCOLS, NROWS.
    add_bounds : bool, default True
        Whether to add cell boundary coordinates.
    chunks : int or dict, optional
        Chunk sizes for the resulting dask-backed dataset.

    Returns
    -------
    xr.Dataset
        The grid dataset.
    """
    gdtyp = metadata["GDTYP"]
    p_alp = metadata["P_ALP"]
    p_bet = metadata["P_BET"]
    xcent = metadata["XCENT"]
    ycent = metadata["YCENT"]
    xorig = metadata["XORIG"]
    yorig = metadata["YORIG"]
    xcell = metadata["XCELL"]
    ycell = metadata["YCELL"]
    ncols = metadata["NCOLS"]
    nrows = metadata["NROWS"]

    if gdtyp == 1:  # Lat-Lon
        crs = "EPSG:4326"
    elif gdtyp == 2:  # Lambert Conformal
        crs = (
            f"+proj=lcc +lat_1={p_alp} +lat_2={p_bet} +lat_0={ycent} "
            f"+lon_0={xcent} +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
        )
    elif gdtyp == 3:  # Mercator
        crs = (
            f"+proj=merc +lat_ts={p_alp} +lon_0={xcent} +lat_0={ycent} "
            f"+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
        )
    elif gdtyp == 4:  # Stereographic
        crs = (
            f"+proj=stere +lat_ts={p_alp} +lat_0={ycent} +lon_0={xcent} "
            f"+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
        )
    elif gdtyp == 5:  # UTM
        crs = f"+proj=utm +zone={int(p_alp)} +datum=WGS84 +units=m +no_defs"
    elif gdtyp == 6:  # Polar Stereographic
        # lat_0 determined by p_alp (1.0 for North, -1.0 for South)
        lat_0 = 90.0 if p_alp > 0 else -90.0
        crs = (
            f"+proj=stere +lat_0={lat_0} +lat_ts={p_bet} +lon_0={xcent} "
            f"+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
        )
    elif gdtyp == 7:  # Equatorial Mercator
        crs = (
            f"+proj=merc +lat_ts={p_alp} +lon_0={xcent} +lat_0=0 "
            f"+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
        )
    elif gdtyp == 8:  # Transverse Mercator
        crs = (
            f"+proj=tmerc +lat_0={ycent} +k={p_bet} +lon_0={xcent} "
            f"+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
        )
    elif gdtyp == 9:  # Albers Equal Area
        crs = (
            f"+proj=aea +lat_1={p_alp} +lat_2={p_bet} +lat_0={ycent} "
            f"+lon_0={xcent} +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
        )
    elif gdtyp == 10:  # Lambert Azimuthal Equal Area
        crs = (
            f"+proj=laea +lat_0={ycent} +lon_0={xcent} "
            f"+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
        )
    elif gdtyp == 13:  # Sinusoidal
        crs = f"+proj=sinu +lon_0={xcent} +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
    else:
        raise ValueError(f"Unsupported IOAPI GDTYP: {gdtyp}")

    extent = (xorig, xorig + ncols * xcell, yorig, yorig + nrows * ycell)
    res = (xcell, ycell)

    ds = create_grid_from_crs(crs, extent, res, add_bounds=add_bounds, chunks=chunks)

    # Attach IOAPI metadata for provenance
    for k, v in metadata.items():
        ds.attrs[f"ioapi_{k}"] = v

    update_history(ds, f"Created grid from IOAPI metadata (GDTYP={gdtyp})")

    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],
    lon_0: float = 0.0,
    radius: float = 6371007.181,
    add_bounds: bool = True,
    chunks: Optional[Union[int, Dict[str, int]]] = None,
) -> xr.Dataset:
    """
    Create a structured grid dataset with a Sinusoidal projection.

    Commonly used for satellite data like MODIS.

    Parameters
    ----------
    extent : Tuple[float, float, float, float]
        Grid extent in projection units (meters): (min_x, max_x, min_y, max_y).
        For global MODIS, this is roughly (-20015109.354, 20015109.354, -10007554.677, 10007554.677).
    res : float or Tuple[float, float] or str
        Grid resolution in meters. If float, same resolution in x and y.
        Supports aliases: '1km' (926.625m), '500m' (463.313m), '250m' (231.656m),
        '5km' (4633.127m), '10km' (9266.254m).
    lon_0 : float, default 0.0
        The central meridian.
    radius : float, default 6371007.181
        Radius of the sphere (default is MODIS sinusoidal sphere radius).
    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'.
    """
    # Resolution aliases for standard MODIS sinusoidal grids
    res_map = {
        "10km": 9266.254331,
        "5km": 4633.1271655,
        "1km": 926.6254331,
        "500m": 463.3127166,
        "250m": 231.6563583,
    }

    if isinstance(res, str) and res in res_map:
        res = res_map[res]

    crs = f"+proj=sinu +lon_0={lon_0} +x_0=0 +y_0=0 +R={radius} +units=m +no_defs"
    ds = create_grid_from_crs(crs, extent, res, add_bounds=add_bounds, chunks=chunks)
    update_history(
        ds, f"Created Sinusoidal grid (lon_0={lon_0}, R={radius}) using xregrid."
    )
    return ds