camsloader.py

"""
Downloads data for a location from CAMS data. Needs cdsapi set-up.
"""
from pathlib import Path
import xarray as xr
import cdsapi
import datetime as dt
import logging
import sys


# Create a logger object with a unique name
logger = logging.getLogger(__name__)
# Set the logging level to INFO (or any other level you prefer)
logger.setLevel(logging.INFO)
# Create a handler to write log messages to stdout
handler = logging.StreamHandler(sys.stdout)
# Set the logging level for the handler
handler.setLevel(logging.INFO)
# Create a formatter to specify the format of the log messages
formatter = logging.Formatter(
    "%(asctime)s - %(name)s - %(levelname)s - %(message)s"
)
# Add the formatter to the handler
handler.setFormatter(formatter)
# Add the handler to the logger object
logger.addHandler(handler)


def dload_cams_data(
    box: list, date: dt.datetime, cams_file: Path | str
) -> tuple:
    """Downloads CAMS data for the region of interest defined in `box`
     and the date in `date`. Units are useful for 6S: AOT550 is
     unitless, TCWV is g/cm2 and O3 is cm-atm.

    Args:
         box (iter): List or similar: max_lat, min_lon, min_lat, max_lon
         date (dt.datetime): date

    Returns:
        tuple: aot550, tcwv, gtco3
    """

    if not cams_file.exists():
        logger.info("Downloading CAMS data for region")
        c = cdsapi.Client()

        c.retrieve(
            "cams-global-reanalysis-eac4",
            {
                "variable": [
                    "black_carbon_aerosol_optical_depth_550nm",
                    "dust_aerosol_optical_depth_550nm",
                    "sea_salt_aerosol_optical_depth_550nm",
                    "sulphate_aerosol_optical_depth_550nm",
                    "total_aerosol_optical_depth_550nm",
                    "total_column_ozone",
                    "total_column_water_vapour",
                ],
                "date": f'{date.strftime("%Y-%m-%d")}/{date.strftime("%Y-%m-%d")}',
                "time": [
                    "00:00",
                    "03:00",
                    "06:00",
                    "09:00",
                    "12:00",
                    "15:00",
                    "18:00",
                    "21:00",
                ],
                "area": box,  # max_lat, min_lon, min_lat, max_lon
                "format": "netcdf",
            },
            cams_file,
        )
    else:
        logger.info("Using cached CAMS data for region")
    f = xr.open_dataset(cams_file, mask_and_scale=True, decode_cf=True)
    f = f.interp(time=date)

    aot550 = f["aod550"][:]  # Unitless
    tcwv = f["tcwv"][:] / 10.0  # g cm^{-2}
    gtco3 = f["gtco3"][:] * 466.98  # Units: cm-atm
    return aot550, tcwv, gtco3