In this paper, a high‐order compact finite difference algorithm is established for the stream function‐velocity formulation of the two‐dimensional steady incompressible Navier‐Stokes equations in general curvilinear coordinates. Different from the previous work, not only the stream function and its first‐order partial derivatives but also the second‐order mixed partial derivative is treated as unknown variable in this work. Numerical examples, including a test problem with an analytical solution, three types of lid‐driven cavity flow problems with unusual shapes and steady flow past a circular cylinder as well as an elliptic cylinder with angle of attack, are solved numerically by the newly proposed scheme. For two types of the lid‐driven trapezoidal cavity flow, we provide the detailed data using the fine grid sizes, which can be considered the benchmark solutions. The results obtained prove that the present numerical method has the ability to solve the incompressible flow for complex geometry in engineering applications, especially by using a nonorthogonal coordinate transformation, with high accuracy.@en