What we believe to be a new rigorous theoretical approach to the refraction of light at the interface of two-dimensional photonic crystals is developed. The proposed method is based on the Dirichlet-to-Neumann (DtN) approach which consists of computing exactly the DtN operators associated with each half-space on both sides of the interface. It fully uses the properties of periodic optical media and takes naturally into account both the evanescent and propagative Bloch modes. Contrary to other proposed approaches, the new method is not based on modal expansions and their complicated electromagnetic field matching at the interfaces, but uses an operator vision. Intrinsically, each operator represents the effect along the interface of a particular medium independently of any medium and/or material that is placed in the other half-space. At the end, the whole computational effort to estimate DtN operators is restricted to the computation of a finite element problem in the periodicity cell of the photonic crystal. Field computations in arbitrary large part of the optical media can be then performed with a negligible computational effort. The method has been applied to the case of incoming plane waves as well as Gaussian beam profiles. It has successfully been compared with the standard plane wave expansion method and finite difference time domain (FDTD) simulations in the case of negative refraction, strongly dispersive, and lensing configurations. The proposed approach is amenable to the generalized study of dispersive phenomena in planar photonic crystals by a rigorous modeling approach avoiding the main drawbacks of FDTD. It is amenable to the study of arbitrary cascaded periodic optical media and photonic crystal heterostructures.