def generate_control_file(fn: str | Path, user_dict: dict, *, simulation_dict: dict[str, Any] | None=None) -> ControlConfig:
    """Generate and write a Cycles control file.

    Provide either direct values or callables that accept a simulation row and return a value in `user_dict`. The
    parameter names should be in lowercase and correspond to the fields in Cycles simulation control files. If a field
    is not provided, it will be filled with a default value. If a field's value is a callable, it will be called with
    the `simulation_dict` to resolve its value.

    The following fields are required in `user_dict`:

      - `simulation_start_year`
      - `simulation_end_year`
      - `rotation_size`
      - `operation_file`
      - `soil_file`
      - `weather_file`

    The default values for other fields are:

      - `crop_file`: `GenericCrops.crop`
      - `reinit_file`: `N/A`
      - `soil_layers`: inferred from the soil file (if not provided)
      - `co2_level`: `-999`
      - `use_reinitialization`: `0`
      - `adjusted_yields`: `0`
      - `hydrology_option`: `1`
      - `automatic_nitrogen`: `0`
      - `automatic_phosphorus`: `0`
      - `automatic_sulfur`: `0`

    All output control fields default to `0`.

    Args:
        fn: Destination control file path.
        user_dict: Values or callables for control fields.
        simulation_dict: Optional simulation row for callable resolution.

    Returns:
        The generated control configuration.
    """
    fn = Path(fn)
    config = _build_control_config(user_dict, simulation_dict, fn.parent)
    write_file(fn, config)

    return config

def read_control_file(control: str | Path) -> ControlConfig:
    """Parse a Cycles control file into a ControlConfig instance.

    Args:
        control: Path to a Cycles control file path.

    Returns:
        Control configuration.
    """
    with open(Path(control)) as f:
        lines = f.read().splitlines()
    lines = iter([line for line in lines if (not line.strip().startswith('#')) and line.strip()])

    hints = get_type_hints(ControlConfig)   # resolves all string annotations → actual types

    control_dict = {}
    for f in fields(ControlConfig):
        target_class = unwrap_optional(hints[f.name])
        sub_hints = get_type_hints(target_class)
        control_dict[f.name] = target_class(**{sub_field.name: parse_value(next(lines), sub_field.name, sub_hints[sub_field.name]) for sub_field in fields(target_class)})

    return ControlConfig(**control_dict)

def generate_nudge_file(fn: str | Path, user_dict: dict, *, simulation_dict: dict[str, Any] | None=None) -> None:
    """Write a Cycles nudge file from user-provided values.

    Provide either direct values or callables that accept a simulation row and return a value in `user_dict`. The
    parameter names should be in lowercase and correspond to the fields in Cycles nudge (calibration) files. If a field
    is not provided, it will be filled with a default value. If a field's value is a callable, it will be called with
    the `simulation_dict` to resolve its value.

    The default values for all calibration multipliers are `1.0`, and the default values for `kd_no3` and `kd_nh4` are
    `0.0` and `5.6`, respectively.

    Args:
        fn: Destination nudge file path.
        user_dict: Values or callables for nudge parameters.
        simulation_dict: Optional simulation row for callable resolution.
    """
    fn = Path(fn)
    config = _build_nudge_config(user_dict, simulation_dict)
    write_file(fn, config)

def read_operation_file(operation: str | Path) -> list:
    """Parse a Cycles operation file into operation objects.

    Args:
        operation: Path to a Cycles operation file.

    Returns:
        A list of operation dataclass instances.
    """
    with open(Path(operation)) as f:
        lines = f.read().splitlines()

    lines = iter([line for line in lines if not line.strip().startswith('#') and line.strip()])
    operations = []
    while True:
        try:
            operation = next(lines).lower()
            if operation not in OPERATION_PARAMETERS:
                raise ValueError(f"Unknown operation keyword found: {operation}")

            target_class = OPERATION_PARAMETERS[operation]
            hints = get_type_hints(target_class)

            operation_dict: dict[str, Any] = {}
            for f in fields(target_class):
                if not f.metadata.get('readable', True):
                    continue
                raw = next(lines)
                if f.name == 'doy':
                    operation_dict[f.name] = raw    # keep as raw string for + detection
                else:
                    operation_dict[f.name] = parse_value(raw, f.name, hints[f.name])

            # Detect relative DOY before constructing the instance
            raw_doy = operation_dict.get('doy', '')
            if isinstance(raw_doy, str) and raw_doy.split()[1].strip().startswith('+'):
                operation_dict['doy'] = int(raw_doy.split()[1].strip()[1:])
                operation_dict['relative_doy'] = True
                operation_dict['resolved_doy'] = _resolve_reference_doy(operations, operation_dict['doy'])
            else:
                operation_dict['doy'] = parse_value(raw_doy, 'doy', hints['doy'])
                operation_dict['relative_doy'] = False

            # Reclassify tillage operations
            if operation == 'tillage':
                tool = operation_dict['tool'].lower().replace('_', '')
                if tool in ('grainharvest', 'harvestgrain'):
                    operation_dict['tool'] = 'grain_harvest'
                    if operation_dict['crop_name'].lower() in ('n/a', 'na', 'all'):
                        operation_dict['crop_name'] = 'All'
                    target_class = Harvest
                elif tool in ('forageharvest', 'harvestforage'):
                    operation_dict['tool'] = 'forage_harvest'
                    if operation_dict['crop_name'].lower() in ('n/a', 'na', 'all'):
                        operation_dict['crop_name'] = 'All'
                    target_class = Harvest
                elif tool in ('kill', 'killcrop', 'killcrops'):
                    operation_dict['tool'] = 'kill'
                    if operation_dict['crop_name'].lower() in ('n/a', 'na', 'all'):
                        operation_dict['crop_name'] = 'All'
                    target_class = Kill

            operations.append(target_class(**operation_dict))
        except StopIteration:
            break

    return operations

def read_output(path: str | Path, output_type: str) -> tuple[pd.DataFrame, dict[str, str]]:
    """Read one Cycles output file and associated unit strings.

    Args:
        path: Directory containing output CSV files.
        output_type: Output file type, such as ``harvest``.

    Returns:
        Parsed DataFrame and a column-to-unit mapping.
    """
    fn = Path(path) / f'{output_type}.csv'
    text = fn.read_text()
    lines = text.splitlines()

    df = pd.read_csv(io.StringIO(text), comment='#').copy()

    for col in DATE_COLUMNS & set(df.columns):
        df[col] = pd.to_datetime(df[col])

    units = _parse_units(lines, df.columns)
    return df, units

def generate_soil_file(fn: str | Path, profile: list[SoilLayer], *,
    target: list[SoilLayer]=DEFAULT_PROFILE, parameters: list[str]=MAPPABLE_PARAMETERS, soil_depth: float | None=None,
    desc: str='', slope: float | None=None, curve_number: float | None=None, hsg: str='') -> list[SoilLayer]:
    """Generate a Cycles-formatted soil file from a list of soil layers.

    The input profile is first mapped to the target layering scheme, then rendered to Cycles soil-file format and
    written to ``fn``.

    Args:
        fn: Output soil file path.
        profile: Soil profile layers in depth order.
        target: Target layer structure to map onto.
        parameters: Soil parameters to map from measured profile.
        soil_depth: Optional maximum depth (m) to include in mapping.
        desc: Optional description/comment line written at file top.
        slope: Optional slope value written to the header.
        curve_number: Optional explicit curve number value.
        hsg: Optional hydrologic soil group used to infer curve number.

    Returns:
        The mapped soil layers written to the output file.

    Raises:
        ValueError: If both ``curve_number`` and ``hsg`` are provided.
    """
    if curve_number is not None and hsg:
        raise ValueError("Only one of curve_number and hsg can be provided.")
    layers = _map_layers(profile, target, parameters, soil_depth)
    Path(fn).write_text('\n'.join(_render_soil_file(layers, desc, slope, curve_number, hsg)) + '\n')
    return layers

def read_soil_file(fn: str | Path) -> tuple[list[SoilLayer], dict]:
    """Read a Cycles soil file into structured objects.

    Args:
        fn: Soil file path to read.

    Returns:
        A tuple containing:
            - List of parsed ``SoilLayer`` objects.
            - Header metadata dictionary with keys such as
              ``curve_number`` and ``slope``.
    """
    lines = [line for line in Path(fn).read_text().splitlines() if line.strip() and not line.strip().startswith('#')]

    meta, data_lines = _parse_header(lines)
    layers = _parse_layers(data_lines)
    return layers, meta

def read_weather_file(fn: str | Path, *, start_year: int=-9999, end_year: int=9999, subdaily: bool=False) -> pd.DataFrame:
    """Read weather file records and filter by year range.

    Args:
        fn: Weather file path.
        start_year: Inclusive first year to keep.
        end_year: Inclusive last year to keep.
        subdaily: If True, parse hourly weather schema.

    Returns:
        Weather data indexed by datetime.
    """
    columns = SUBDAILY_COLUMNS if subdaily else DAILY_COLUMNS
    df = pd.read_csv(
        Path(fn),
        usecols=list(range(len(columns))),
        names=list(columns.keys()),
        comment='#',
        sep=r'\s+',
        na_values='-999',
        skiprows=WEATHER_HEADER_LINES,
    )
    df = df.astype(columns)
    df.index = _build_date_index(df, subdaily=subdaily)
    df.index.name = 'date'
    return df[(df['YEAR'] >= start_year) & (df['YEAR'] <= end_year)]

def generate_reinit_file(out_path: str | Path, in_path: str | Path, doy: int) -> None:
    """Write a re-initialization file from Cycles output for all years at a selected day-of-year.

    Args:
        out_path: Destination re-initialization file path.
        in_path: Simulation output directory that contains Cycles output files.
        doy: Day-of-year extracted from re-initialization.
    """
    out_path = Path(out_path)
    in_path  = Path(in_path)

    output_df, _ = read_output(in_path, 'reinit')

    n_layers = sum(1 for col in output_df.columns if col.startswith('soil_moisture_content'))

    output_df['year'] = output_df['date'].dt.year
    filtered = output_df[output_df['date'].dt.dayofyear == doy].reset_index(drop=True)

    lines: list[str] = []
    for _, row in filtered.iterrows():
        lines.extend(_format_reinit_block(row, doy, n_layers))

    out_path.write_text('\n'.join(lines) + '\n')

def plot_yield(harvest_df: pd.DataFrame, *, ax: Axes | None=None, fontsize: int | None=None) -> Axes:
    """Plot grain and forage yields by crop.

    Args:
        harvest_df: Harvest output DataFrame.
        ax: Optional axes to draw on.
        fontsize: Optional global font size override.

    Returns:
        Axes containing the yield plot.
    """
    if ax is None:
        _, ax = plt.subplots()
    if fontsize is not None:
        plt.rcParams.update({'font.size': fontsize})

    crops = harvest_df['crop'].unique().tolist()
    crop_colors = _assign_crop_colors(crops, ax)

    for crop in crops:
        for harvest, marker in HARVEST_MARKERS.items():
            _plot_harvest_type(ax, harvest_df, crop, harvest, marker, crop_colors[crop])

    ax.set_ylabel(YIELD_UNIT_LABEL)
    ax.set_axisbelow(True)
    ax.grid(True, color='#93a1a1', alpha=0.2)
    ax.legend(
        handles=_build_legend_handles(crops, crop_colors),
        handletextpad=0,
        bbox_to_anchor= (1.0, 0.5),
        loc='center left',
        shadow=True,
        frameon=False,
    )
    return ax

def plot_operations(operations: list, rotation_size: int, *, axs: Axes | np.ndarray | None=None, fontsize: int | None=None):
    """Plot operations by day-of-year for each rotation year.

    Args:
        operations: Sequence of parsed operation objects.
        rotation_size: Number of years in the rotation.
        axs: Optional axes object(s) for rendering.
        fontsize: Optional global font size override.

    Returns:
        Axes array used to render timelines.
    """
    if axs is None:
        _, axs = plt.subplots(rotation_size, 1, sharex=True)
    assert axs is not None

    if isinstance(axs, Axes):
        axs = np.array(axs).reshape((1,))

    if rotation_size != axs.shape[0]:
        raise ValueError('The number of axes must match the rotation size.')

    if fontsize is not None: plt.rcParams.update({'font.size': fontsize})

    for y in range(rotation_size):
        for key, value in OPERATION_TYPES.items():
            sub_list = [op for op in operations if type(op).__name__.lower() == key and op.year == y + 1]

            if not sub_list: continue

            axs[y].plot(
                [op.doy if not op.relative_doy else op.resolved_doy for op in sub_list], [value.yloc] * len(sub_list),
                'o',
                label=value.title + ':\n' + '\n'.join(f'{op.doy if not op.relative_doy else op.resolved_doy}: {getattr(op, value.label)}' if value.label is not None else f'{op.doy}' for op in sub_list),
                color=value.color,
                ms=10,
            )

            if key == 'planting':
                for op in sub_list:
                    if op.end_doy == -999: continue
                    axs[y].fill_betweenx([value.yloc - 0.5, value.yloc + 0.5], op.doy, op.end_doy, alpha=0.25, color=value.color)

        axs[y].set_xlim(-1, 370)
        axs[y].grid(False)
        axs[y].yaxis.set_ticks_position('none')
        axs[y].yaxis.set_tick_params(left=False, right=False, which='both', labelleft=False)
        axs[y].set_ylim(-3, 7)
        axs[y].text(184, 5, f'Year {y + 1}', ha='center')
        for spine in ('right', 'left', 'top'):
            axs[y].spines[spine].set_color('none')

        # set the y-spine
        axs[y].spines['bottom'].set_position('zero')

        # turn off the top spine/ticks
        axs[y].xaxis.tick_bottom()
        axs[y].set_xticks(MDOYS)
        axs[y].set_xticklabels(MONTHS)

        handles, _ = axs[y].get_legend_handles_labels()
        if handles:
            axs[y].legend(
                loc='center left',
                bbox_to_anchor=(1.1, 0.5),
                ncols=5,
                frameon=False,
            )

    return axs

def plot_map(gdf: gpd.GeoDataFrame, column: str, *, projection: ccrs.Projection=ccrs.PlateCarree(), ax: Sequence[float] | GeoAxes | None=None,
    cmap: Colormap | str='viridis', vmin: float | None=None, vmax: float | None=None,
    colorbar: bool=True, cb_axes: tuple[float, float, float, float] | None=None, extend: str='neither', cb_orientation: str='horizontal',
    label: str | None=None, title: str | None=None,
    fontsize: float | None=None,
    frameon: bool=False) -> tuple[Figure, GeoAxes]:
    """Render a thematic map from a GeoDataFrame column.

    Args:
        gdf: GeoDataFrame to visualize.
        column: Column name to visualize.
        projection: Map projection for output axes.
        ax: Existing GeoAxes or add_axes rectangle.
        cmap: Matplotlib colormap.
        vmin: Optional lower bound for colormap normalization.
        vmax: Optional upper bound for colormap normalization.
        colorbar: Whether to draw a colorbar.
        cb_axes: Optional colorbar axes rectangle.
        extend: Colorbar extension mode.
        cb_orientation: Colorbar orientation.
        label: Optional colorbar label.
        title: Optional plot title.
        fontsize: Optional global font size override.
        frameon: Whether to draw map frame and grid labels.

    Returns:
        Tuple of figure and GeoAxes.
    """

    if fontsize is not None: plt.rcParams.update({'font.size': fontsize})

    if ax is None:
        fig = plt.figure(figsize=(9, 6))
        ax = fig.add_axes((0.025, 0.09, 0.95, 0.93), projection=projection, frameon=frameon)    # type: ignore
    elif isinstance(ax, Sequence):
        fig = plt.figure(figsize=(9, 6))
        ax = fig.add_axes(ax, projection=projection, frameon=frameon)   # type: ignore
    elif isinstance(ax, GeoAxes):
        fig = ax.get_figure()

    if colorbar is True:
        cax = fig.add_axes((0.3, 0.07, 0.4, 0.02) if cb_axes is None else cb_axes)  # type: ignore

    gdf.plot(
        column=column,
        cmap=cmap,
        ax=ax,  # type: ignore
        vmin=vmin,
        vmax=vmax,
    )
    ax.add_feature(feature.STATES, edgecolor=[0.7, 0.7, 0.7], linewidth=0.5)    # type: ignore
    ax.add_feature(feature.LAND, facecolor=[0.8, 0.8, 0.8])     # type: ignore
    ax.add_feature(feature.LAKES)   # type: ignore
    ax.add_feature(feature.OCEAN)   # type: ignore

    if frameon:
        gl = ax.gridlines(      # type: ignore
            draw_labels=True,
            color='gray',
            dms=True,
            x_inline=False,
            y_inline=False,
            linestyle='--',
        )
        gl.bottom_labels = None # type: ignore
        gl.right_labels = None  # type: ignore

    if colorbar is True:
        cbar = plt.colorbar(
            ax.collections[0],  # type: ignore
            cax=cax,
            orientation=cb_orientation,
            extend=extend,
        )
        if label is not None: cbar.set_label(label)
        cbar.ax.xaxis.set_label_position('top' if cb_orientation == 'horizontal' else 'right')  # type: ignore
    if title is not None:
        ax.set_title(title) # type: ignore

    return fig, ax  # type: ignore

def plot_satellite_map(fig: Figure, extent: tuple[float, float, float, float], *,
    alpha: float=1.0, ax: tuple[float, float, float, float] | None=None, desired_pixels: int=1024, style: str | None=None) -> tuple[GeoAxes, ccrs.Projection]:
    """Render a satellite basemap over a geographic extent.

    Args:
        fig: Matplotlib figure used to create the map axes.
        extent: Geographic bounds as ``(west, east, south, north)`` in degrees.
        alpha: Basemap transparency where ``1.0`` is fully opaque.
        ax: Optional axes rectangle passed to ``Figure.add_axes``.  If omitted, a full-figure subplot is created.
        desired_pixels: Target figure resolution in pixels used to estimate tile zoom.
        style: Optional Cartopy GoogleTiles style string. If omitted, ArcGIS World Imagery tiles are used.

    Returns:
        Tuple containing:
            - GeoAxes with the satellite basemap.
            - Projection associated with the tile source.
    """
    GOOGLE_URL = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.jpg'

    google_satellite = cimgt.GoogleTiles(url=GOOGLE_URL) if style is None else cimgt.GoogleTiles(style=style)

    if ax is None:
        ax = fig.add_subplot(1, 1, 1, projection=google_satellite.crs)  # type: ignore
    else:
        ax = fig.add_axes(ax, projection=google_satellite.crs)  # type: ignore

    assert isinstance(ax, GeoAxes)
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    ax.add_image(google_satellite, _zoom_from_extent(extent, desired_pixels=desired_pixels), alpha=alpha)

    return ax, google_satellite.crs