In this work, we consider the two-dimensional Navier--Stokes equations with boundary conditions on the velocity on one part of the boundary and involving the pressure fluid on the rest of the boundary. We write a new weak formulation in terms of the stream function of the velocity field. Next, we propose C1-virtual element methods of high-order k \geq 2 to discretize this formulation, and a novel analysis to prove optimal error estimates in the H2-norm under minimal regularity conditions on the weak stream function is developed. Moreover, algorithms to compute the velocity, pressure, and vorticity fields as a postprocess of the discrete stream function are proposed. Furthermore, we provide a comparison between our methods and other numerical schemes in terms of computational standpoint, which indicates that the stream function virtual element formulation is an attractive, efficient, and competitive alternative to solve the Navier--Stokes system. Finally, we report several benchmark numerical tests to support our theoretical findings and highlight interesting features of our discrete scheme.