The interaction between planetary waves and an arbitrary zonal flow is studied from a phase-space viewpoint. Using the Wigner distribution, a planetary wave Vlasov equation is derived that includes the contribution of the mean flow to the zonal potential vorticity gradient. This equation is applied to the problem of planetary wave modulational instability, where it is used to predict a fastest growing mode of finite wavenumber. A wave-mean flow numerical model is used to test the analytical predictions, and an intuitive explanation of modulational instability and jet asymmetry is given via the motion of planetary wavepackets in phase space.