For linear perturbations of black hole spacetimes, the ringdown is described by quasinormal modes. Quasinormal modes are resonance states, defined by ingoing and outgoing radiation conditions at the horizon and infinity. Their wavefunctions do not lie in a Hilbert space, and they do not in general form a complete basis. In particular, there is a priori no obvious inner product under which they are orthonormal. In contrast to normal modes, this limits efforts to carry out higher order perturbation theory in terms of quasinormal modes. In this work, we show that for type D spacetimes with a t-φ reflection symmetry, gravitational quasinormal modes are in fact orthogonal. We work in terms of Weyl scalars. On this space we define a symmetric bilinear form ⟨⟨·, ·⟩⟩ in terms of the symplectic form for the Teukolsky equation on a Cauchy surface Σ and the t-φ reflection operator. (Our bilinear form is motivated by earlier work of Leung, Liu, and Young [Phys. Rev. A 49, 3057 (1994)].) The bilinear form is complex-linear in both entries, it is independent of the precise choice of Σ, and the time-evolution operator is symmetric with respect to ⟨⟨·, ·⟩⟩. It follows that quasinormal modes are orthogonal with respect to ⟨⟨·, ·⟩⟩ and have finite norm. We also relate the bilinear form on Weyl scalars to bilinear forms on the spaces of metric perturbations and Hertz potentials. The formula for a quasinormal mode excitation coefficient emerges naturally as the projection of initial data onto the quasinormal mode, however this projector also extends to data consisting of quasinormal modes. By projecting with the bilinear form, our goal is to develop a framework for higher order black hole perturbation theory in terms of quasinormal modes.