Galactic Equator zenith transit
If I eliminate ec and z, and do this:

    eq = ephem.Equatorial(ra, dec, epoch = d)
    zg = ephem.Galactic(eq)


the crossings I find in gzt.csv agree with Stellarium within a minute of time.

If you're open to another package by the same developer, skyfield.almanac can help with such calculations.
In this case you would implement a function which takes a time argument and returns an integer indicating which side of the galactic equator the zenith is on.
Then you would pass this function as an argument to find_discrete.