Strong-form meshfree methods for incompressible flows represent boundary shapes using discrete points (DPs) and enable local evaluation of governing equations without global meshes. However, such methods suffer from checkerboard instability due to the collocated arrangement of velocity and pressure variables. Moreover, in projection-based meshfree incompressible flow solvers using node-based moving least-squares (MLS), constructing a discrete Laplacian operator consistent with the divergence of the gradient is difficult, leading to insufficient enforcement of the solenoidal condition. A staggered MLS is proposed in [1] introducing a local primal–dual grid, which enables compatible evaluation of discrete gradient and divergence operators. However, projection-based meshfree incompressible flow solvers that satisfy the incompressibility constraint at the discrete level have not yet been established. We propose a projection-based strong-form meshfree method for incompressible flows based on the staggered MLS [2], built upon a mesh-constrained discrete point (MCD) method [3,4]. In the MCD method, each DP is associated one-to-one with a background mesh element, and is allowed to exist only within its corresponding mesh element (mesh constraint). This mesh constraint enables accurate representation of shapes using DPs and efficient stencil computation through structured mesh-based indexing. For spatial discretization, a staggered MLS based on the local primal–dual grid proposed in [1] is employed. The radial velocity component is evaluated on interfaces of virtual dual cells and used for MLS reconstruction. In addition, a time evolution converting formula [5] is employed to temporally link velocity and pressure variables and enable strong coupling. Numerical tests confirm that the proposed method strictly satisfies the solenoidal condition of the velocity field at the discrete level, that the flow field is more stably calculated than the collocated method at higher time resolution, and that the proposed method can represent the flow characteristics more accurately than the collocated method even at higher Reynolds numbers.