I believe the answer is positive.

First step: reduce to the case that $\mathrm{rank}(df_x)\geq n-1$ except on a finite union of submanifolds of dimension at most $n-4.$

*Proof:* For $r=0,\dots,n-2,$ let $M_r\subset \mathbb R^{n\times n}$ be the set of matrices of rank exactly $r.$ Each of these is a smooth manifold of dimension at most $n^2-4.$ For each such $r,$ consider the map $\phi_r:\mathrm{int}(\mathbb D^n)\times M_r\to\mathbb R^{n\times n}$ defined by $\phi_r(x,M)=M-df_x.$ A preimage $\phi_r^{-1}(\{N\})$ consists of pairs $(x,df_x+N)$ where $x$ is such that $df_x+N$ has rank $r$ i.e. such that $x\mapsto f(x)+Nx$ has a derivative of rank $r.$
By the preimage theorem, if $N$ is regular for $\phi_r$ then this preimage is a manifold of codimension at least $4.$ Since this preimage is actually a graph of the smooth function $x\mapsto df_x+N,$ projecting away the second component doesn't affect the codimension. A similar argument can be applied to the boundary of $\mathbb D^n.$ By Sard's theorem, almost all $N$ are regular values for all these maps. Take a small such $N$ so that $x\mapsto f(x)+Nx$ is approximately $f.$ $\square$

Actually it will be useful to also require $\mathrm{rank}(df_0)\geq n-1,$ which is easy enough because it holds for almost all $N$ - just pick $N$ outside $(\bigcup_r M_r)-df_0.$

Second step: given $\epsilon>0,$ we can construct a diffeo $\phi$ from $\mathrm{int}(\mathbb D^n)$ to a subset of the points $x$ with $\mathrm{rank}(df_x)\geq n-1,$ such that $\max_x|d\phi_x|= O(1/\epsilon),$ and the set where $\phi(x)\neq x$ has measure $O(\epsilon^3).$ The implicit constants here are allowed to depend on $f.$

*Proof*: Let $E$ denote the exceptional set $\{x\in\mathrm{int}(\mathbb D^n)\mid \mathrm{rank}(df_x)<n-1\}.$ The idea is to push $E$ away from the origin in a radial direction, out of the disc.

Let $\epsilon>0.$ Define $$R_0=\min_{x\in E}|x|.$$ Define $E'$ to be the set of points $v\in S^{n-1}$ such that $\lambda v\in E$ for some $0\leq \lambda < 1$ - the radial projection of $E$ to the boundary of the disc.
I will argue that the $\epsilon$-neighborhood of $E'$ in $S^{n-1}$ (the set of points at distance at most $\epsilon$ from $E'$) has $(n-1)$-dimensional volume at most $O(\epsilon^3).$
The homeomorphism $\phi:(v,r)\mapsto rv$ from $S^{n-1}\times (R_0/2,1)$ to $\{x\in D^n\mid |x|\in(R_0/2,1)\}$ is a diffeo. So the preimage of the the at-most-$(n-4)$-dimensional manifold $E$ still has *Minkowski* dimension $n-4.$
It can therefore be covered by a set of $N(\epsilon)$ balls $B((v_i,r_i),\epsilon)$ such that $N(\epsilon) \epsilon^{n-4}=O(1).$ The projected balls $B(v_i,\epsilon)\subset S^{n-1}$ cover $E'$ and still satisfy $N(\epsilon) \epsilon^{n-4}=O(1).$ This means the $(n-1)$-dimensional volume is at most $O(1)N(\epsilon) \epsilon^{n-1}=O(\epsilon^3).$
Doubling the radii gives a set $\bigcup_i B(v_i,2\epsilon)$ still of $(n-1)$-dimensional volume $O(\epsilon^3),$ and containing the $\epsilon$-neighborhood of $E'.$

We will need a function $h:S^{n-1}\to[R_0/4,2]$ such that:

- $h$ is smooth
- $h\leq R_0/2$ on $E'$
- $h\geq 1$ except on the $\epsilon$-neighborhood of $E'$
- $h$ has Lipschitz constant $O(1/\epsilon)$

To construct such a function, first define $g(v)=R_0/4$ for $v\in E',$ and $g(v)=2$ for $v$ not in the $\epsilon$-neighborhood of $E'.$
Note that $g$ has Lipschitz constant at most $2/\epsilon.$
By the Kirszbraun theorem, $g$ extends to a function $S^{n-1}\to [R_0/4,2]$ with the same Lipschitz constant.

However $g$ may not be smooth. To deal with this, we can think of $S^{n-1}$ as the symmetric space $SO(n)/SO(n-1),$ so $g$ is a function on $SO(n)$ invariant under the right action of $SO(n-1)$ (specifically, $SO(n-1)$ is the stabilizer of $(1,0,0,\dots,0),$ and a rotation $\rho\in SO(n)$ corresponds to the point $\rho\cdot (1,0,0,\dots,0)\in S^{n-1}$). Convolve by left multiplication by a smooth mollifier $SO(n)\to\mathbb R$ supported on a small ball around the identity. The resulting function is still invariant under the right $SO(n-1)$ action, so corresponds to a function $h:S^{n-1}\to [R_0/4,2].$ This function $h$ will be smooth (basically by the implicit function theorem). It is a convex combination of rotated versions of $g,$ where each rotation moves points by no more than $R_0\epsilon/8$ say. This ensures conditions (2.) and (3.), and also (4.) because the set of functions with Lipschitz constant $\leq 2/\epsilon$ is convex.

Let $\psi:\mathbb R\to\mathbb R$ be a smooth function with $\psi(x)=x$ for $x\leq 0,$ and $\psi(x)<R_0/2$ and $0<\psi'(x)\leq 1$ everywhere.
Define $\phi:\mathrm{int}(\mathbb D^n)\to \mathrm{int}(\mathbb D^n)$ by $\phi(0)=0$ and $$\phi(rv)=(\psi(r-h(v))+h(v))v$$ for unit vectors $v$ and real $0<r<1.$
Then $\phi(rv)=rv$ if $r\leq h(v),$ in particular for $r\leq R_0/4$ or $h(v)\geq 1.$
So $\phi(x)=x$ in a neighborhood of zero, which ensures $\phi$ is smooth at zero. And $\phi(rv)=rv$ for all $0<r<1$ and $v\not\in E'$; since $E'$ has $(n-1)$-dimensional measure $O(\epsilon^3),$ the set of $x$ of the form $rv$ with $v\not\in E'$ has $n$-dimensional measure $O(\epsilon^3).$ (The probability that a uniformly chosen $x\in D^n$ satisfies $x/|x|\in E'$ is the same as the probability that a uniformly chosen $v\in S^{n-1}$ lies in $E'.$)

The Lipschitz constant of $\phi$ can be bounded by bounding the derivatives of $\phi$ along a path through $(r,v)$ with $r>R_0/4$ and $|\dot r|,|\dot v|\leq 1.$
Writing $h=h(v)$ and $\psi=\psi(r-h),$ the chain rule gives $|\dot h|=O(1/\epsilon)$ and hence $|\dot\psi|=O(1/\epsilon).$
So $$\frac{d}{dt} \phi(rv) = (\dot\psi+\dot h)v + (\psi+h)\dot v = O(1/\epsilon)$$
which means $\max_x|d\phi_x|=O(1/\epsilon).$ $\square$

The integral of $|d(f\circ \phi)-df|^2$ is at most the maximum value of the integrand, $O(1/\epsilon^2),$ multiplied by the volume of the set where the integrand can be non-zero, $O(\epsilon^3),$ giving $O(\epsilon)$ overall. The integral of $|(f\circ \phi)-f|^2$ is even better: $O(\epsilon^3).$ So $\|f-f\circ \phi\|^2_{W^{1,2}}=O(\epsilon).$ And the rank of $f\circ \phi$ will be $\geq n-1$ everywhere.