API Reference

def __init_subclass__(cls, database_path: str, **kwargs):
    """Automatically sets the database path for subclasses.

    Parameters
    ----------
    database_path : str or None
        The directory path where the specific atmosphere data is stored.
        Pass ``None`` for sources that do not use a file database
        (e.g. ``CustomAtmosphere`` or API-based retrievers).
    """
    super().__init_subclass__(**kwargs)
    cls.database_path = None if database_path is None else Path(database_path)

@classmethod
def at_sites(
    cls,
    times: np.ndarray[tuple[int], np.dtype[np.datetime64]] | pd.DatetimeIndex,
    latitude: Sequence[float] | float,
    longitude: Sequence[float] | float,
    site_names: Sequence[str] | None = None,
) -> Self:
    """Load atmospheric data at specific geographic locations.

    Extracts and interpolates MERRA-2 data for one or more point locations.
    Performs bilinear spatial interpolation and quadratic temporal interpolation.

    Parameters
    ----------
    times : np.ndarray or pd.DatetimeIndex
        Time points for data extraction.
    latitude : float or Sequence[float]
        Latitude coordinate(s) in decimal degrees.
    longitude : float or Sequence[float]
        Longitude coordinate(s) in decimal degrees.
    site_names : Sequence[str] or None, default None
        Optional names for the sites.

    Returns
    -------
    MERRA2DailyAtmosphere
        Instance containing interpolated atmospheric data.

    Raises
    ------
    ValueError
        If latitude and longitude have different lengths or invalid values.

    Examples
    --------
        >>> import pandas as pd
        >>> from spartasolar.atmoslib import MERRA2DailyAtmosphere

        >>> # Single location
        >>> times = pd.date_range("2020-01-01", periods=5, freq="D")
        >>> atmos = MERRA2DailyAtmosphere.at_sites(
        ...     times=times,
        ...     latitude=40.4168,
        ...     longitude=-3.7038,
        ...     site_names="Madrid"
        ... )

        >>> # Multiple locations
        >>> lats = [40.4168, 41.3851, 36.7213]  # Madrid, Barcelona, Málaga
        >>> lons = [-3.7038, 2.1734, -4.4214]
        >>> names = ["Madrid", "Barcelona", "Málaga"]
        >>> atmos = MERRA2DailyAtmosphere.at_sites(
        ...     times=times,
        ...     latitude=lats,
        ...     longitude=lons,
        ...     site_names=names
        ... )

        >>> # Access data
        >>> pressure = atmos.dataset["pressure"]
        >>> print(pressure.dims)
        ('time', 'site')
    """

    latitude = [latitude] if isinstance(latitude, (float, int)) else latitude
    latitude = np.asarray([validate_type(lat, Latitude) for lat in latitude], dtype=float).reshape(-1)
    longitude = [longitude] if isinstance(longitude, (float, int)) else longitude
    longitude = np.asarray([validate_type(lon, Longitude) for lon in longitude], dtype=float).reshape(-1)

    if len(latitude) != len(longitude):
        raise ValueError('latitude and longitude must have the same length')

    # load the dataset. Check for local availability. If not available, download.
    dataset = cls._load_dataset(times)

    # lat-lon interpolation
    output_lat = xr.DataArray(latitude, dims="site", name="lat")
    output_lon = xr.DataArray(longitude, dims="site", name="lon")
    # # if interpolation was `nearest` it is faster to use .sel than .interp
    # output_dataset = dataset.sel(lat=output_lat, lon=output_lon, method='nearest')
    # output_dataset = dataset.interp(lat=output_lat, lon=output_lon, method='nearest')
    output_dataset = dataset.interp(lat=output_lat, lon=output_lon, method='linear')

    # time interpolation
    if "time" in output_dataset.coords:
        output_dataset = (
            output_dataset
            .assign_coords(
                time=pd.to_datetime(output_dataset.time.values).tz_localize(None)  # convert to tz-naive for interpolation
            )
            .interp(
                time=ensure_tz_aware_datetime_index(times, utc=True).tz_localize(None),  # convert to tz-naive for interpolation
                method='quadratic'
            )
        )

    global_attrs = {
        "title": "Daily Clear-sky Atmospheric Dataset for SPARTA",
        "references": "doi:10.5067/KLICLTZ8EM9D, doi:10.5067/Q9QMY5PBNV1T, doi:10.5067/VJAFPLI1CSIV",
    }

    obj = cls()
    obj._atmosphere = build_atmosphere_of_sites(
        times=times,
        latitude=latitude,
        longitude=longitude,
        constituents=output_dataset.data_vars,
        site_names=site_names,
        global_attrs=global_attrs)
    return obj

def compute(
    self,
    model: Model = "SPARTA",
    include_atmosphere: bool = False,  # TODO: currently ignored
    model_kwargs: dict | None = None,
) -> xr.Dataset:
    """Compute clear-sky solar radiation using a radiative transfer model.

    This method integrates solar position calculations with atmospheric
    constituent data to compute clear-sky irradiance components (GHI, DNI,
    DHI, etc.) using the specified radiative transfer model.

    Parameters
    ----------
    model : Model, default "SPARTA"
        Name of the clear-sky model to use. Options: "SPARTA", "Bird"
    include_atmosphere : bool, default False
        If True, include atmospheric constituents in the output dataset.
        If False, only radiation components are returned.
    model_kwargs : dict, optional
        Additional keyword arguments to pass to the model function

    Returns
    -------
    xr.Dataset
        CF-compliant dataset containing computed irradiance components:
        - ghi: Global Horizontal Irradiance (W/m²)
        - dni: Direct Normal Irradiance (W/m²)
        - dhi or dif: Diffuse Horizontal Irradiance (W/m²)
        - csi: Circumsolar Irradiance (W/m², SPARTA only)

    Examples
    --------
    >>> import pandas as pd
    >>> from spartasolar.atmoslib import MERRA2DailyAtmosphere
    >>>
    >>> times = pd.date_range("2020-06-15", periods=24, freq="h")
    >>> atm = MERRA2DailyAtmosphere.at_sites(
    ...     times=times,
    ...     latitude=36.72,
    ...     longitude=-4.42
    ... )
    >>> result = atm.compute(model="SPARTA")
    >>> print(result.ghi.values)

    Use different model with custom parameters:

    >>> result = atm.compute(
    ...     model="Bird",
    ...     model_kwargs={"scheme": "transmittance_parameterization"}
    ... )

    Notes
    -----
    The method automatically:
    - Calculates solar position (zenith angle, Earth-Sun distance)
    - Converts atmospheric units to model requirements
    - Handles both gridded and site-based data structures

    See Also
    --------
    spartasolar.modlib.sparta : SPARTA model implementation
    spartasolar.modlib.bird : Bird clear-sky model
    """

    model = validate_type(model, Model)
    model_func = getattr(modlib, model)
    model_vars = inspect.getfullargspec(model_func).args

    is_regular_grid = "site" not in self.dataset.dims

    # compute solar geometry...
    sw_eval = sunwhere.regular_grid if is_regular_grid else sunwhere.sites
    solpos = sw_eval(
        # self.dataset.time.values,  # WARNING: this returns a utc tz-naive numpy array of dtype datetime64[ns] !!!
        self.dataset.indexes["time"],  # NOTE: with `indexes` we get a pandas DatetimeIndex with tz info preserved!
        latitude=self.dataset.lat.values,
        longitude=self.dataset.lon.values,
        algorithm=config.get_option("sunwhere.algorithm", default="psa"),
        refraction=config.get_option("sunwhere.refraction", default=True),
        engine=config.get_option("sunwhere.engine", default="numexpr")
    )

    new_coord_names = {"location": "site", "latitude": "lat", "longitude": "lon"}
    cosz = solpos.cosz.rename({old: new for old, new in new_coord_names.items() if old in solpos.cosz.dims})

    if is_regular_grid:
        n_lats = self.dataset.sizes["lat"]
        n_lons = self.dataset.sizes["lon"]
        ecf = solpos.ecf.expand_dims(dim={"lat": n_lats, "lon": n_lons}, axis=(1, 2))
    else:
        n_locs = self.dataset.sizes["site"]
        ecf = solpos.ecf.expand_dims(dim={"site": n_locs}, axis=1)

    kwargs = {"cosz": cosz, "ecf": ecf}

    # atmosphere...
    def get_with_proper_units(var):
        if var not in self.dataset.data_vars:
            raise ValueError(f"variable `{var}` required by the model `{model}` is not available in this atmosphere")
        if var == "pressure":
            return self.dataset[var] * 1e-2  # Pa to hPa
        if var == "pwater":
            return pwater_in_kg_m2_to_cm(self.dataset[var])
        if var == "ozone":
            return ozone_in_kg_m2_to_cm(self.dataset[var])
        return self.dataset[var]
    variables = set(model_vars).intersection(self.dataset.data_vars)  # variables required by this model
    kwargs = kwargs | {var: get_with_proper_units(var) for var in variables}

    # and call the clearsky model...
    result = model_func(**(kwargs | (model_kwargs or {})))

    # encapsulate the result in a CF-compliant xarray Dataset
    if is_regular_grid:
        return build_atmosphere_on_regular_grid(
            times=self.dataset.indexes["time"],
            latitude=self.dataset.lat.values,
            longitude=self.dataset.lon.values,
            constituents=result,
            global_attrs=get_global_attrs(feature_type="grid")
        )
    else:
        site_values = self.dataset.coords.get("site").values
        n_sites = self.dataset.sizes["site"]
        site_names = validate_site_names(site_values, n_sites)
        return build_atmosphere_of_sites(
            times=self.dataset.indexes["time"],
            latitude=self.dataset.lat.values,
            longitude=self.dataset.lon.values,
            constituents=result,
            site_names=site_names,
            global_attrs=get_global_attrs(feature_type="timeSeries")
        )

@classmethod
def on_regular_grid(
    cls,
    times: np.ndarray[tuple[int], np.dtype[np.datetime64]] | pd.DatetimeIndex,
    latitude: Sequence[float] | float,
    longitude: Sequence[float] | float,
) -> Self:
    """Load atmospheric data on a regular lat/lon grid.

    Extracts and interpolates MERRA-2 data onto a user-specified regular grid.
    The output has dimensions (time, lat, lon). NaN values in albedo are filled
    with zeros.

    Parameters
    ----------
    times : np.ndarray or pd.DatetimeIndex
        Time points for data extraction.
    latitude : float or Sequence[float]
        Latitude coordinates for the grid in decimal degrees.
    longitude : float or Sequence[float]
        Longitude coordinates for the grid in decimal degrees.

    Returns
    -------
    MERRA2DailyAtmosphere
        Instance containing gridded atmospheric data.

    Raises
    ------
    ValueError
        If coordinates are out of valid range.

    Examples
    --------
        >>> import pandas as pd
        >>> import numpy as np
        >>> from spartasolar.atmoslib import MERRA2DailyAtmosphere

        >>> # Define a regular grid over Iberian Peninsula
        >>> times = pd.date_range("2020-06-01", periods=10, freq="D")
        >>> lats = np.linspace(36, 44, 9)  # 9 latitude points
        >>> lons = np.linspace(-10, 4, 15)  # 15 longitude points

        >>> atmos = MERRA2DailyAtmosphere.on_regular_grid(
        ...     times=times,
        ...     latitude=lats,
        ...     longitude=lons
        ... )

        >>> # Access gridded data
        >>> albedo = atmos.dataset["albedo"]
        >>> print(albedo.dims)
        ('time', 'lat', 'lon')
        >>> print(albedo.shape)
        (10, 9, 15)
    """

    latitude = [latitude] if isinstance(latitude, (float, int)) else latitude
    latitude = np.asarray([validate_type(lat, Latitude) for lat in latitude], dtype=float).reshape(-1)
    longitude = [longitude] if isinstance(longitude, (float, int)) else longitude
    longitude = np.asarray([validate_type(lon, Longitude) for lon in longitude], dtype=float).reshape(-1)

    logger.info(f"loading MERRA-2 daily data for {len(times)} time points, "
                f"{len(latitude)} latitudes, {len(longitude)} longitudes...")

    # load the dataset. Check for local availability. If not available, download.
    dataset = cls._load_dataset(times)
    dataset["albedo"] = dataset["albedo"].fillna(0.)

    # lat-lon interpolation
    output_dataset = dataset.interp(lat=latitude, lon=longitude, method='linear')

    # time interpolation
    if "time" in output_dataset.coords:
        output_dataset = (
            output_dataset
            .assign_coords(
                time=pd.to_datetime(output_dataset.time.values).tz_localize(None)  # convert to tz-naive for interpolation
            )
            .interp(
                time=ensure_tz_aware_datetime_index(times, utc=True).tz_localize(None),  # convert to tz-naive for interpolation
                method='quadratic'
            )
        )

    global_attrs = {
        "title": "Daily Clear-sky Atmospheric Dataset for SPARTA",
        "references": "doi:10.5067/KLICLTZ8EM9D, doi:10.5067/Q9QMY5PBNV1T, doi:10.5067/VJAFPLI1CSIV",
    }

    obj = cls()
    obj._atmosphere = build_atmosphere_on_regular_grid(
        times=times,
        latitude=latitude,
        longitude=longitude,
        constituents=output_dataset.data_vars,
        global_attrs=global_attrs)
    return obj

def __init_subclass__(cls, database_path: str, **kwargs):
    """Automatically sets the database path for subclasses.

    Parameters
    ----------
    database_path : str or None
        The directory path where the specific atmosphere data is stored.
        Pass ``None`` for sources that do not use a file database
        (e.g. ``CustomAtmosphere`` or API-based retrievers).
    """
    super().__init_subclass__(**kwargs)
    cls.database_path = None if database_path is None else Path(database_path)

@classmethod
def at_sites(
    cls,
    times: np.ndarray[tuple[int], np.dtype[np.datetime64]] | pd.DatetimeIndex,
    latitude: Sequence[float] | float,
    longitude: Sequence[float] | float,
    site_names: Sequence[str] | None = None,
) -> Self:
    """Retrieve monthly climatology at specific sites.

    Parameters
    ----------
    times : np.ndarray or pd.DatetimeIndex
        Time stamps for climatology retrieval. Monthly climatology is
        repeated for each year and interpolated to exact times.
    latitude : Sequence[float]
        Latitude(s) in degrees North [-90, 90]
    longitude : Sequence[float]
        Longitude(s) in degrees East [-180, 180]
    site_names : Sequence[str], optional
        Names for each site

    Returns
    -------
    MERRA2LTAAtmosphere
        Instance with interpolated climatological data

    Examples
    --------
    >>> times = pd.date_range("2023-01-01", "2023-12-31", freq="D")
    >>> atm = MERRA2LTAAtmosphere.at_sites(
    ...     times=times,
    ...     latitude=36.72,
    ...     longitude=-4.42,
    ...     site_names="Málaga"
    ... )
    """

    latitude = np.asarray(latitude, dtype=float).reshape(-1)
    latitude = [validate_type(lat, Latitude) for lat in latitude]
    longitude = np.asarray(longitude, dtype=float).reshape(-1)
    longitude = [validate_type(lon, Longitude) for lon in longitude]

    if len(latitude) != len(longitude):
        raise ValueError('latitude and longitude must have the same length')

    # load the dataset. Check for local availability. If not available, download.
    dataset = cls._load_dataset(times)

    # lat-lon interpolation
    output_lat = xr.DataArray(latitude, dims="site", name="lat")
    output_lon = xr.DataArray(longitude, dims="site", name="lon")
    output_dataset = dataset.interp(lat=output_lat, lon=output_lon, method='linear')

    # time interpolation
    if "time" in output_dataset.coords:
        output_dataset = (
            output_dataset
            .assign_coords(
                time=pd.to_datetime(output_dataset.time.values).tz_localize(None)  # convert to tz-naive for interpolation
            )
            .interp(
                time=ensure_tz_aware_datetime_index(times, utc=True).tz_localize(None),  # convert to tz-naive for interpolation
                method='quadratic'
            )
        )

    global_attrs = {
        "title": "Long Term Average Clear-sky Atmospheric Dataset for SPARTA",
        "references": "doi:10.5067/KLICLTZ8EM9D, doi:10.5067/Q9QMY5PBNV1T, doi:10.5067/VJAFPLI1CSIV",
    }

    obj = cls()
    obj._atmosphere = build_atmosphere_of_sites(
        times=times,
        latitude=latitude,
        longitude=longitude,
        constituents=output_dataset.data_vars,
        site_names=site_names,
        global_attrs=global_attrs)
    return obj

def compute(
    self,
    model: Model = "SPARTA",
    include_atmosphere: bool = False,  # TODO: currently ignored
    model_kwargs: dict | None = None,
) -> xr.Dataset:
    """Compute clear-sky solar radiation using a radiative transfer model.

    This method integrates solar position calculations with atmospheric
    constituent data to compute clear-sky irradiance components (GHI, DNI,
    DHI, etc.) using the specified radiative transfer model.

    Parameters
    ----------
    model : Model, default "SPARTA"
        Name of the clear-sky model to use. Options: "SPARTA", "Bird"
    include_atmosphere : bool, default False
        If True, include atmospheric constituents in the output dataset.
        If False, only radiation components are returned.
    model_kwargs : dict, optional
        Additional keyword arguments to pass to the model function

    Returns
    -------
    xr.Dataset
        CF-compliant dataset containing computed irradiance components:
        - ghi: Global Horizontal Irradiance (W/m²)
        - dni: Direct Normal Irradiance (W/m²)
        - dhi or dif: Diffuse Horizontal Irradiance (W/m²)
        - csi: Circumsolar Irradiance (W/m², SPARTA only)

    Examples
    --------
    >>> import pandas as pd
    >>> from spartasolar.atmoslib import MERRA2DailyAtmosphere
    >>>
    >>> times = pd.date_range("2020-06-15", periods=24, freq="h")
    >>> atm = MERRA2DailyAtmosphere.at_sites(
    ...     times=times,
    ...     latitude=36.72,
    ...     longitude=-4.42
    ... )
    >>> result = atm.compute(model="SPARTA")
    >>> print(result.ghi.values)

    Use different model with custom parameters:

    >>> result = atm.compute(
    ...     model="Bird",
    ...     model_kwargs={"scheme": "transmittance_parameterization"}
    ... )

    Notes
    -----
    The method automatically:
    - Calculates solar position (zenith angle, Earth-Sun distance)
    - Converts atmospheric units to model requirements
    - Handles both gridded and site-based data structures

    See Also
    --------
    spartasolar.modlib.sparta : SPARTA model implementation
    spartasolar.modlib.bird : Bird clear-sky model
    """

    model = validate_type(model, Model)
    model_func = getattr(modlib, model)
    model_vars = inspect.getfullargspec(model_func).args

    is_regular_grid = "site" not in self.dataset.dims

    # compute solar geometry...
    sw_eval = sunwhere.regular_grid if is_regular_grid else sunwhere.sites
    solpos = sw_eval(
        # self.dataset.time.values,  # WARNING: this returns a utc tz-naive numpy array of dtype datetime64[ns] !!!
        self.dataset.indexes["time"],  # NOTE: with `indexes` we get a pandas DatetimeIndex with tz info preserved!
        latitude=self.dataset.lat.values,
        longitude=self.dataset.lon.values,
        algorithm=config.get_option("sunwhere.algorithm", default="psa"),
        refraction=config.get_option("sunwhere.refraction", default=True),
        engine=config.get_option("sunwhere.engine", default="numexpr")
    )

    new_coord_names = {"location": "site", "latitude": "lat", "longitude": "lon"}
    cosz = solpos.cosz.rename({old: new for old, new in new_coord_names.items() if old in solpos.cosz.dims})

    if is_regular_grid:
        n_lats = self.dataset.sizes["lat"]
        n_lons = self.dataset.sizes["lon"]
        ecf = solpos.ecf.expand_dims(dim={"lat": n_lats, "lon": n_lons}, axis=(1, 2))
    else:
        n_locs = self.dataset.sizes["site"]
        ecf = solpos.ecf.expand_dims(dim={"site": n_locs}, axis=1)

    kwargs = {"cosz": cosz, "ecf": ecf}

    # atmosphere...
    def get_with_proper_units(var):
        if var not in self.dataset.data_vars:
            raise ValueError(f"variable `{var}` required by the model `{model}` is not available in this atmosphere")
        if var == "pressure":
            return self.dataset[var] * 1e-2  # Pa to hPa
        if var == "pwater":
            return pwater_in_kg_m2_to_cm(self.dataset[var])
        if var == "ozone":
            return ozone_in_kg_m2_to_cm(self.dataset[var])
        return self.dataset[var]
    variables = set(model_vars).intersection(self.dataset.data_vars)  # variables required by this model
    kwargs = kwargs | {var: get_with_proper_units(var) for var in variables}

    # and call the clearsky model...
    result = model_func(**(kwargs | (model_kwargs or {})))

    # encapsulate the result in a CF-compliant xarray Dataset
    if is_regular_grid:
        return build_atmosphere_on_regular_grid(
            times=self.dataset.indexes["time"],
            latitude=self.dataset.lat.values,
            longitude=self.dataset.lon.values,
            constituents=result,
            global_attrs=get_global_attrs(feature_type="grid")
        )
    else:
        site_values = self.dataset.coords.get("site").values
        n_sites = self.dataset.sizes["site"]
        site_names = validate_site_names(site_values, n_sites)
        return build_atmosphere_of_sites(
            times=self.dataset.indexes["time"],
            latitude=self.dataset.lat.values,
            longitude=self.dataset.lon.values,
            constituents=result,
            site_names=site_names,
            global_attrs=get_global_attrs(feature_type="timeSeries")
        )

@classmethod
def on_regular_grid(
    cls,
    times: np.ndarray[tuple[int], np.dtype[np.datetime64]] | pd.DatetimeIndex,
    latitude: Sequence[float],
    longitude: Sequence[float],
) -> Self:
    """Retrieve monthly climatology on a regular spatial grid.

    Parameters
    ----------
    times : np.ndarray or pd.DatetimeIndex
        Time stamps for climatology retrieval
    latitude : Sequence[float]
        Latitude grid coordinates in degrees North
    longitude : Sequence[float]
        Longitude grid coordinates in degrees East

    Returns
    -------
    MERRA2LTAAtmosphere
        Instance with gridded climatological data

    Examples
    --------
    >>> import numpy as np
    >>> lats = np.linspace(36.0, 41.0, 20)
    >>> lons = np.linspace(-5.0, -3.0, 20)
    >>> times = pd.date_range("2023-01-15", periods=12, freq="MS") + pd.Timedelta(14.5, "D")
    >>> atm = MERRA2LTAAtmosphere.on_regular_grid(
    ...     times=times,
    ...     latitude=lats,
    ...     longitude=lons
    ... )
    """

    latitude = np.asarray(latitude, dtype=float).reshape(-1)
    latitude = [validate_type(lat, Latitude) for lat in latitude]
    longitude = np.asarray(longitude, dtype=float).reshape(-1)
    longitude = [validate_type(lon, Longitude) for lon in longitude]

    # load the dataset. Check for local availability. If not available, download.
    dataset = cls._load_dataset(times)

    # lat-lon interpolation
    output_dataset = dataset.interp(lat=latitude, lon=longitude, method='linear')

    # time interpolation
    if "time" in output_dataset.coords:
        output_dataset = (
            output_dataset
            .assign_coords(
                time=pd.to_datetime(output_dataset.time.values).tz_localize(None)  # convert to tz-naive for interpolation
            )
            .interp(
                time=ensure_tz_aware_datetime_index(times, utc=True).tz_localize(None),  # convert to tz-naive for interpolation
                method='quadratic'
            )
        )

    global_attrs = {
        "title": "Long Term Average Clear-sky Atmospheric Dataset for SPARTA",
        "references": "doi:10.5067/KLICLTZ8EM9D, doi:10.5067/Q9QMY5PBNV1T, doi:10.5067/VJAFPLI1CSIV",
    }

    obj = cls()
    obj._atmosphere = build_atmosphere_on_regular_grid(
        times=times,
        latitude=latitude,
        longitude=longitude,
        constituents=output_dataset.data_vars,
        global_attrs=global_attrs)
    return obj

def __init_subclass__(cls, database_path: str, **kwargs):
    """Automatically sets the database path for subclasses.

    Parameters
    ----------
    database_path : str or None
        The directory path where the specific atmosphere data is stored.
        Pass ``None`` for sources that do not use a file database
        (e.g. ``CustomAtmosphere`` or API-based retrievers).
    """
    super().__init_subclass__(**kwargs)
    cls.database_path = None if database_path is None else Path(database_path)

@classmethod
def at_site(
    cls,
    times: pd.DatetimeIndex | np.ndarray[tuple[int], np.dtype[np.datetime64]],
    latitude: float,
    longitude: float,
    site_name: str | None = None,
) -> Self:
    """Retrieve MERRA-2 data from GEE for a specific site.

    Downloads data from GEE if not cached, applies spatial/temporal
    corrections, then interpolates to requested times.

    Parameters
    ----------
    times : pd.DatetimeIndex or np.ndarray of datetime64
        Time stamps for data retrieval (UTC)
    latitude : float
        Latitude in degrees North [-90, 90]
    longitude : float
        Longitude in degrees East [-180, 180]
    site_name : str, optional
        Name identifier for the site

    Returns
    -------
    MERRA2GEEAtmosphere
        Instance with interpolated atmospheric data

    Examples
    --------
    >>> import pandas as pd
    >>> from spartasolar import config
    >>> config.set_option('merra2_gee.project', 'my-gee-project')
    >>>
    >>> times = pd.date_range("2020-06-15", periods=24, freq="h")
    >>> atm = MERRA2GEEAtmosphere.at_site(
    ...     times=times,
    ...     latitude=36.72,
    ...     longitude=-4.42,
    ...     site_name="Málaga"
    ... )
    >>> result = atm.compute(model="SPARTA")

    Notes
    -----
    - Requires `merra2_gee.project` configuration
    - Automatically corrects GEE's 0.25° latitude offset
    - Adjusts time stamps from hour-start to hour-center
    - Data is cached locally to minimize API calls
    """

    latitude = validate_type(latitude, Latitude)
    longitude = validate_type(longitude, Longitude)

    def fetch_and_distill_and_archive(year: int, path: Path) -> None:
        logger.info(f"fetching GEE data for year={int(year)} and path={path.as_posix()}")
        if not (gee_project := get_option("merra2_gee.project")):
            raise ValueError(
                'missing Google cloud\'s project. Add `project = "<your_gee_project>"` '
                f"in the `merra2_gee` table in `{get_config_path()}` and reload spartasolar "
                "or use `spartasolar.config.set_option('merra2_gee.project', <your_project>)`"
            )
        logger.info(f"using GEE project <green>{gee_project}</green>")
        data = fetch_merra2_data_from_gee_api(
            latitude=latitude,
            longitude=longitude,
            start_date=f"{year}-01-01",
            end_date=f"{year}-12-31",
            project=gee_project,
        )

        data = cls._distill_crude_data(data, latitude, longitude)
        logger.debug(f"{data.head()=}")
        data.to_parquet(path)
        logger.success(f"data downloaded and archived in <blue>{path.name}</blue>")

    # load data from one year before and one year after the requested times_utc, but
    # clipping the years on 2004 and the current year
    paths = []
    times_utc = ensure_tz_aware_datetime_index(times, utc=True)
    years = cls._infer_years_from_times(times_utc)
    for year in sorted(years):
        if not (path := cls._get_filename(year, latitude, longitude)).exists():
            fetch_and_distill_and_archive(year, path)
        paths.append(path)
    logger.debug([path.as_posix() for path in paths])
    data = pd.read_parquet(paths)  # WARNING: times in data are UTC, tz-naive !!!

    # interpolate to times_utc
    x = times_utc.to_numpy(dtype="datetime64[s]").astype(float)
    xi = data.times_utc.to_numpy(dtype="datetime64[s]").astype(float)
    data_i = data.drop(columns=["times_utc", "albedo"])
    y = interp1d(xi, data_i.values, kind="quadratic", axis=0)(x)
    y_dict = {key: value for key, value in zip(data_i.columns, y.T)}
    y_dict["albedo"] = interp1d(xi, data.albedo.values, kind="linear", axis=0)(x)
    data_interp = pd.DataFrame({"times": times_utc} | y_dict)

    global_attrs = {
        "title:": "HourlyEarth Engine Data Catalog dataset for SPARTA",
        "source": "NASA/GMAO MERRA-2 reanalysis via Google Earth Engine API",
        "references": "doi:10.5067/KLICLTZ8EM9D, doi:10.5067/Q9QMY5PBNV1T, doi:10.5067/VJAFPLI1CSIV",
    }

    obj = cls()
    obj._atmosphere = build_atmosphere_of_sites(
        times=times,  # WARNING: return the input times, even if they are not tz-aware!
        latitude=latitude,
        longitude=longitude,
        constituents=data_interp.drop(columns=["times"]),
        site_names=site_name,
        global_attrs=global_attrs,
    )
    return obj