I'm using GDAL python to create a bunch of tiff files from a netcdf file. 

    def modelout_to_geotif(netcdfile, outdir, VarName):
        """ Given the path to the netcdf output,
        produces geotiffs for each band of the given variable name.
        
        Input::
            netcdfile: string path to netcdf file
            VarName: string variable name in netcdf file (e.g. 'GridPrec')
            
        Output::
            creates geotiffs for each 'band' (e.g. time step). Output is written
            to current working directory.
            
        
        """
        try:
            nci = gdal.Open('NETCDF:{0}:{1}'.format(netcdfile, VarName))
            ncd = Dataset(netcdfile)
        except:
            print "Could not open input file: {0}".format(netcdfile)
            return
        
        try:
            geotransform = nci.GetGeoTransform()
        except:
            print "Could not get geotransform.\n"
            return
        try:
            projection = nci.GetProjection()
        except:
            print "Could not get projection.\n"
            return
        
        numbands = nci.RasterCount
        times = [o2d(t) for t in ncd.variables['time']]
        print "Found {0} bands to process for {1}".format(numbands, VarName)
        for i in range(numbands):
            time = times[i]
            band = nci.GetRasterBand(i+1)
            raster = band.ReadAsArray()
            y,x = raster.shape
            
            output_file = "{0}_{1}.tif".format(VarName, time.strftime('%Y%m%d_%H'))
            subdir = os.path.join(outdir, VarName)
            if not os.path.exists(subdir):
                os.mkdir(subdir)
            output_file = os.path.join(subdir, output_file)
            if os.path.exists(output_file):
                continue
            # Create gtif
            driver = gdal.GetDriverByName("GTiff")
            #TODO: Make this more generic
            dst_ds = driver.Create(output_file, x, y, 1, gdal.GDT_CFloat32 )
            
            # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
            dst_ds.SetGeoTransform( geotransform )
      
            # set the reference info 
            dst_ds.SetProjection( projection )
            
            # write the band
            dst_ds.GetRasterBand(1).WriteArray(raster)


Now I would like to add those to a mosaic catalog in an ESRI File Geodatabase. I've managed to get it working using arcpy, but it's really slow, and I'd prefer to use GDAL. Is there a way? Below is my arcpy script:

    import os
    try:
        import arcpy
        from arcpy import env
    except:
        raise ImportError("Unable to import arcpy - are you using ESRI Python?")
    
    env = 'P:/Enki'
    
    
    def add_to_mosaic(filegdb, raster_dir):    
        
        print env
        #Add Raster Dataset type Raster to FGDB Mosaic Dataset
        #Calculate Cell Size Ranges and Build Boundary
        #Build Overviews for Mosaic Dataset upon the 3rd level Raster Dataset pyramid
        #Apply TIFF file filter
        #Build Pyramids for the source datasets
        varnames = ['GridPrec']#,
                    #'GridSnowOut',
                    #'GridRad',
                    #'GridTemp',
                    #'GridSoilOut',
                    #'GridRunoff']
        
        for VarName in varnames:
            gdbname = filegdb 
            mdname = os.path.join(gdbname, VarName)
            messages = 'Failed'
            sr =  arcpy.SpatialReference()    
    #        sr.name = 'WGS_1984_UTM_Zone_32N'
            sr.factoryCode = 32632
            sr.create()
            print sr 
            try:
                noband = 1
                res = None
                #TODO: get from input data
                print 'trying to create %s' % VarName
                res = arcpy.CreateMosaicDataset_management(gdbname, mdname, sr, noband)

            except:
                print res.getMessages()
                
                
            
            
            rastype = "Raster Dataset"
            inpath = os.path.join(raster_dir, VarName)
            updatecs = "UPDATE_CELL_SIZES"
            updatebnd = "UPDATE_BOUNDARY"
            updateovr = "UPDATE_OVERVIEWS"
            maxlevel = "2"
            maxcs = "#"
            maxdim = "#"
            spatialref = "#"
            inputdatafilter = "*200701*.tif"
            subfolder = "NO_SUBFOLDERS"
            duplicate = "EXCLUDE_DUPLICATES"
            buildpy = "BUILD_PYRAMIDS"
            calcstats = "CALCULATE_STATISTICS"
            buildthumb = "NO_THUMBNAILS"
            comments = "Add Raster Datasets"
            forcesr = "#"
        
            res = arcpy.AddRastersToMosaicDataset_management(mdname,  rastype, inpath, updatecs, \
                                                       updatebnd, updateovr, maxlevel, maxcs, \
                                                       maxdim, spatialref, inputdatafilter, \
                                                       subfolder, duplicate, buildpy, calcstats, \
                                                       buildthumb, comments, forcesr)
            print res.getMessages()