We assume that we observe N independent copies of a diffusion process on a time-interval [0,2T]. For a given time t, we estimate the transition density p_t(x,y), namely the conditional density of X_t + s given X_s = x, under conditions on the diffusion coefficients ensuring that this quantity exists. We use a least squares projection method on a product of finite dimensional spaces, prove risk bounds for the estimator and propose an anisotropic model selection method, relying on several reference norms. A simulation study illustrates the theoretical part for Ornstein-Uhlenbeck or square-root (Cox-Ingersoll-Ross) processes.