def load_training_rasters(response_raster, explanatory_rasters, selected=None):
    """
    Parameters
    ----------
    response_raster : Path to GDAL raster containing responses
    explanatory_rasters : List of Paths to GDAL rasters containing explanatory variables

    Returns
    -------
    train_xs : Array of explanatory variables
    train_ys : 1xN array of known responses
    """

    with rasterio.open(response_raster) as src:
        response_data = src.read().flatten()

    if selected is None:
        train_y = response_data
    else:
        train_y = response_data[selected]

    selected_data = []
    for rast in explanatory_rasters:
        with rasterio.open(rast) as src:
            explanatory_data = src.read().flatten()
        assert explanatory_data.size == response_data.size
        if selected is None:
            selected_data.append(explanatory_data)
        else:
            selected_data.append(explanatory_data[selected])

    train_xs = np.asarray(selected_data).T
    return train_xs, train_y