I like the discussion (only possible in this dimension) which uses (1) the fact that Calderón—Zygmund operators which are smoothing of order one transform bounded measurable functions into continuous ones with $e\log1/e$ modulus of continuity and (2) Osgood's elementary theorem that this modulus of continuity is good enough for unique solutions of ODEs.
Now in a chart with a bounded measurable change of almost complex structure thought of at each point as point in the unit disc with non-Euclidean geometry and at the origin the standard structure in the chart, draw the Poincaré geodesic between these two structures.
Applying the Cauchy transform converts the bounded measurable infinitesimal distortion of conformal structure into an Osgood vector field.
We can integrate this and move along the geodesic path.
This produces a composition of homeomorphisms with the required bounded conformal distortion.
To get charts with holomorphic overlap mappings modify the starting charts by these constructed homeomorphisms with bounded conformal distortion and use this fact: "homeomorphisms with bounded conformal distortion and zero distortion ae. are holomorphic".
This proves Newlander—Nirenberg for bounded measurable almost complex structures (relative to a standard background).
This is not the most elementary proof of the smooth result but to me it is the conceptually easiest proof.
It also gives a natural and very strong statement with lots of applications not possible using the smooth result.
Also the same proof scheme (use Calderón—Zygmund then Osgood to inch your way to a solution) also solves the Euler equation for 2D incompressible fluid motion for any fixed time which is not known in higher dimensions.