Conventional mathematical models for simulating incompressible fluid flow problems are based on the Navier-Stokes equations expressed in terms of pressure and velocity. In this context, pressure-velocity coupling is a key issue, and countless numerical techniques and methods have been developed over the decades to solve these equations efficiently and accurately. In two dimensions, an alternative approach is to rewrite the Navier-Stokes equations regarding two scalar quantities: the streamfunction and the vorticity. Compared to the primitive variables approach, this formulation does not require pressure to be computed, thereby avoiding the inherent difficulties associated with the pressure-velocity coupling. However, deriving boundary conditions for the streamfunction and vorticity is challenging. This work proposes an efficient, high-order accurate finite-volume discretisation of the two-dimensional incompressible Navier-Stokes equations in the streamfunction-vorticity formulation. A detailed discussion is devoted to deriving the appropriate boundary conditions and their numerical treatment, including on arbitrary curved boundaries. The reconstruction for off-site data method is employed to avoid the difficulties associated with generating curved meshes to preserve high-orders of convergence in arbitrary curved domains, such as sophisticated meshing algorithms, cumbersome quadrature rules, and intricate non-linear transformations. This method approximates arbitrary curved boundaries with a conventional linear piecewise approximation, while constrained polynomial reconstructions near the boundary fulfil the prescribed conditions at the physical boundary. Several incompressible fluid flow test cases in non-trivial 2D curved domains are presented and discussed to demonstrate the accuracy and effectiveness of the proposed methodology in achieving very high orders of convergence.