We consider upscaling of non-Darcy flow in heterogeneous porous media. Our approach extends the pressure-based numerical homogenization procedure for linear Darcy flow, due to Durlofsky, to the nonlinear case. The effective coefficients are not constants but rather mildly varying functions of prevailing gradients of pressure. The upscaled model approximates the fine grid model accurately and, in some cases, more accurately than what is expected for Darcy flow; this is due to the non-Darcy effects which suppress heterogeneity. We provide comparisons of alternative approaches as well as consider several variants of numerical realizations of the non-Darcy flow model. Numerical results show effectiveness of the upscaling procedure.