A weak Galerkin meshless (WGM) method is proposed and analyzed in this paper for the incompressible stationary Navier–Stokes equations. In Galerkin meshless methods, we have to deal with the integration of non-polynomial functions, and then the often-used high-order Gauss quadrature rules not only lead to high computational cost but also seriously damage the optimal convergence. In the WGM method, a weakly defined gradient is introduced to facilitate numerical integration and restore the optimal convergence. Stability of the WGM method is analyzed, optimal order error estimates of the velocity and the pressure are derived, and convergence of the Oseen iteration for dealing with the nonlinear convection term is proved. Theoretical results reveal the basic principle of selecting quadrature rules in meshless methods to ensure that the optimal convergence is completely independent of numerical integration. Numerical results show the performance of the proposed WGM method and verify theoretical results.