In this paper we address the large-scale regularity theory for the stationary Navier-Stokes equations in highly oscillating bumpy John domains. These domains are very rough, possibly with fractals or cusps, at the microscopic scale, but are amenable to the mathematical analysis of the Navier-Stokes equations. We prove: (i) a large-scale Calderon-Zygmund estimate, (ii) a large-scale Lipschitz estimate, (iii) large-scale higher-order regularity estimates, namely, $C^{1,gamma}$ and $C^{2,gamma}$ estimates. These nice regularity results are inherited only at mesoscopic scales, and clearly fail in general at the microscopic scales. We emphasize that the large-scale $C^{1,gamma}$ regularity is obtained by using first-order boundary layers constructed via a new argument. The large-scale $C^{2,gamma}$ regularity relies on the construction of second-order boundary layers, which allows for certain boundary data with linear growth at spatial infinity. To the best of our knowledge, our work is the first to carry out such an analysis. In the wake of many works in quantitative homogenization, our results strongly advocate in favor of considering the boundary regularity of the solutions to fluid equations as a multiscale problem, with improved regularity at or above a certain scale.