Assuming the CSV XY values form a regular grid, you can treat the CSV as an [XYZ raster][1], 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)



  [1]: https://gdal.org/drivers/raster/xyz.html