How to perform spatial join and obtain attribute in python/GDAL?

I am working on a script that will take an unprojected raster, and using a UTM grid (shapefile), will obtain the correct SRS (I have a shapefile where each zone is attributed with its numerical SRS) and perform a warp/projection using GDAL on that raster. It is a quicker way than figuring out the right zone, and then manually doing the warp.

Procedurally, I am grabbing the raster corners and constructing WKT geometry (polygon); the objective is to do a spatial join with the correct feature in the UTM grid. Here's what I have so far:

import os

import ogr
from osgeo import gdal

def getRasterBox(raster):
    src = gdal.Open(raster)
    ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform()
    lrx = ulx + (src.RasterXSize * xres)
    lry = uly + (src.RasterYSize * yres)
    return (uls, uly, lrx, lry)


def findUTMzone(raster_box, utmgrid):
    # do the stuff

utmzone_path = r'C:\temp'
utmzone_file = 'utm_grid.shp'
utmzone = os.path.join(utmzone_path, utmzone_file)

utm_ds = ogr.Open(utmzone)
utm_layer = utm_ds.GetLayer()


rasterdir = r'C:\temp'
raster_file = 'n38.tif'
raster = os.path.join(rasterdir, raster_file)

rasterdetails = getRasterBox(raster)
rasterWKT = 'POLYGON(({} {}, {} {}, {} {}, {} {}, {} {}))'.format(rasterdetails[0], rasterdetails[1], rasterdetails[0], rasterdetails[3], rasterdetails[2], rasterdetails[3], rasterdetails[2], rasterdetails[1], rasterdetails[0], rasterdetails[1])

Using the QuickWKT plugin for QGIS, I can see that I'm getting a good WKT ring for the raster bounds. However, I'm unsure how to perform the spatial join in a way that grabs the feature attribute containing the correct SRS code. How do I do this part (under def findUTMzone(raster_box, utmgrid):)? I have seen a few arcpy examples, but am looking for something open-source.