The density distribution of hard spheres in square and rectangular channels is studied using density functional theory and grand canonical ensemble Monte Carlo simulations. The theory uses the weighted density approximation for the excess Helmholtz free energy functional with the Curtin-Ashcroft approximation for the weighting function. The theoretical predictions are in good quantitative agreement with the simulations except for very high densities. The theory predicts pronounced layering in the channel as the density is increased, with high values of density at the surfaces and even higher densities at the corners where any two surfaces meet. Interesting interference effects are observed in the density profiles when compared to the case of a hard sphere fluid in slitlike pores.