import gdal
import ogr

stream = ogr.Open("stream.shp")
substrates = ogr.Open("substrate1.shp")
stream_l = stream.GetLayer()
sub_l = substrates.GetLayer()

driver = ogr.GetDriverByName("ESRI Shapefile")
outputds = driver.CreateDataSource(outputfn)
outputlayer = outputds.CreateLayer(outputfn, geom_type=ogr.wkbPolygon)

for stream_f  in stream_l:
    s_geom = stream_f.GetGeometryRef()

   #Get overlap size by using a clip
    dict = {}
    for i,sub_feat in enumerate(sub_l):
        geom = sub_feat.GetGeometryRef()
        clip = geom.Intersection(stream_geom)
        dict[clip.GetArea()]=i

    #Sort based on size and get most overlapping polygon
    sorted(dict.items(), key=lambda s: s[0])
    polynum = dict[dict.keys[0]]
    largest_feat = sub_l.GetFeature(polynum)

    #enter the fieldname at #fn
    field_largest_feat = largest_feat.GetField(#fn)

    #add this substrate feature to the original stream feature
    defn = outputlayer.GetLayerDefn()
    feature =ogr.Feature(defn)
    feature.CreateField(ogr.FieldDefn("substrate",ogr.OFTString))
    feature.SetField(str(field_largest_feat))
    feature.SetGeometry(s_geom)

    outputlayer.CreateFeature(feature)

outputds = None