The Koopman operator approach provides a powerful linear description of nonlinear dynamical systems in terms of the evolution of observables. While the operator is typically infinite-dimensional, it is crucial to develop finite-dimensional approximation methods and characterize the related approximation errors with upper bounds, preferably expressed in the uniform norm. In this paper, we depart from the traditional use of orthogonal projection or truncation, and propose a novel method based on Bernstein polynomial approximation. Considering a basis of Bernstein polynomials, we construct a matrix approximation of the Koopman operator in a computationally effective way. Building on results of approximation theory, we characterize the rates of convergence and the upper bounds of the error in various contexts including the cases of univariate and multivariate systems, and continuous and differentiable observables. The obtained bounds are expressed in the uniform norm in terms of the modulus of continuity of the observables. Finally, the method is extended to a data-driven setting through a proper change of coordinates. Numerical experiments show that the method is robust to noise and demonstrates good performance for trajectory prediction.