Assuming the CSV XY values form a regular grid, you can treat the CSV as an XYZ raster, though you will need to extract out just the X, Y and specific "Z" column from the CSV for each band that will form the multiband output. You can then clip to your shapefile with gdal.Warp and stack the clipped results.

import pandas as pd
from osgeo import gdal
from urllib.request import urlretrieve
from zipfile import ZipFile

# # Download file from ESRI
# url = (
#        'https://www.arcgis.com/sharing/rest/content/items/715db3ed501b42fe9581caaa5c56caf9/data'
#        )
# filename = 'countyshapefiles.zip'
# urlretrieve(url, filename) # the zip has a size of ~80 mb
#
# # Unzip
# with ZipFile('countyshapefiles.zip', 'r') as z:
#     z.extractall(
#         path = 'countyshapefiles'
#     )

shapefile = r"countyshapefiles/USA_Counties.shp"
epsg = 32617
csvfile = "hill_multicolumn.csv"  # Assumes X, Y, col1,...,colN

df = pd.read_csv(csvfile)
df = df.sort_values(["Y", "X"], ascending=True)
vrts = []
for col in df.columns[2:]:
    df_col = df.filter(["X", "Y", col])
    df_col.to_csv(f"hill_{col}.csv", index=False)
    rds = gdal.OpenEx(f"hill_{col}.csv", gdal.OF_RASTER)
    vrt = gdal.Warp(
        f"/vsimem/hill_{col}.vrt",
        rds,
        cutlineDSName=shapefile,
        cutlineWhere="NAME = 'Hillsborough' AND STATE_NAME = 'Florida'",
        cropToCutline=True,
        srcSRS=f"EPSG:{epsg}"
    )
    vrts.append(vrt)

vrt = gdal.BuildVRT(f"/vsimem/hill.vrt", vrts)
multiband = gdal.Translate("hill_multiband.tif", vrt)

Assuming the CSV XY values form a regular grid, you can treat the CSV as an XYZ raster, though you will need to extract out just the X, Y and specific "Z" column from the CSV for each band that will form the multiband output. You can then clip to your shapefile with gdal.Warp and stack the clipped results.

import pandas as pd
from osgeo import gdal
from urllib.request import urlretrieve
from zipfile import ZipFile

# # Download file from ESRI
# url = (
#        'https://www.arcgis.com/sharing/rest/content/items/715db3ed501b42fe9581caaa5c56caf9/data'
#        )
# filename = 'countyshapefiles.zip'
# urlretrieve(url, filename) # the zip has a size of ~80 mb
#
# # Unzip
# with ZipFile('countyshapefiles.zip', 'r') as z:
#     z.extractall(
#         path = 'countyshapefiles'
#     )

shapefile = r"countyshapefiles/USA_Counties.shp"
epsg = 32617
csvfile = "hill_multicolumn.csv"  # Assumes X, Y, col1,...,colN

df = pd.read_csv(csvfile)
df = df.sort_values(["Y", "X"], ascending=True)
vrts = []
for col in df.columns[2:]:
    df_col = df.filter(["X", "Y", col])
    df_col.to_csv(f"hill_{col}.csv", index=False)
    rds = gdal.OpenEx(f"hill_{col}.csv", gdal.OF_RASTER)
    vrt = gdal.Warp(
        f"/vsimem/hill_{col}.vrt",
        rds,
        cutlineDSName=shapefile,
        cutlineWhere="NAME = 'Hillsborough' AND STATE_NAME = 'Florida'",
        cropToCutline=True,
        srcSRS=f"EPSG:{epsg}"
    )
    vrts.append(vrt)

vrt = gdal.BuildVRT(f"/vsimem/hill.vrt", vrts, separate=True)
multiband = gdal.Translate("hill_multiband.tif", vrt)