Prescribing a gradient direction Let $\Omega= \{(x,y) : \frac{1}{2} \leq x^2+y^2 \leq 1\}$ and $S = \{(x,y) :  x^2+y^2 = 1\}$ the unit circle, and $X=w^{1.\infty}(\Omega;\mathbb{R})$ the space of Lipschitz valued functions. We denote by $\left| \cdot \right|$ the euclidean norm in $\mathbb{R}^2$ or the modulus in $\mathbb{R}$.
Given $\phi \in C^\infty(\Omega;S)$, I am looking for a continuous linear operator $P$ acting on $X$ such that 
For all $f\in X$,
$$
|D(Pf)\cdot \phi| = |D(Pf)| ,
$$
and
$$\text{ If } |Df \cdot \phi| = |Df| \mbox{ in } \Omega, \text{ then }P(f)=f.$$
In other words, if the gradient of  $f$ is parallel to $\phi$ then $f$ is unchanged by $P$, whereas if it isn't, $P$ "projects" it on that space. 
There is one simple case when one can do it naturally, it is when $\phi = (x,y)/|(x,y)|$, as in this case one can choose
$$
Pf(x,y) = \frac{1}{2\pi}\int_{0}^{2\pi} f\left(|(x,y)|\cos \theta, |(x,y)| \sin \theta \right) d\theta.
$$
As $Pf$ is radial,   $D(Pf)$ is parallel to $\phi$.  
If there exist two sufficiently smooth functions such that $u,v$ such that 
a. $\phi= Du$ 
b. $\min \det (Du,Dv)>0$ 
c. $v$ is such that for all $u$
$$
\int_{\{u=c\}\cap\{x\in \Omega\}} dv=1
$$
then, writing 
$$
x=s(u,v) \text{ and } y=t(u,v)
$$
we can generalize the radial case to
$$
Pf(x,y) = \int_{\{v : (s(u,v),t(u,v))\in \Omega \}} f\left(s(u(x,y),v),t(u(x,y),v)\right) dv.
$$
As pointed out by  Mikhail Skopenkov in his remark, this is easier said than done. Suppose $\phi=(1,0)$, then $u=x$; but the natural choice of $v=y/l(x)$, where
$$
l(x)=\begin{cases} 2\sqrt{1-x^2} & \text{ when } \frac{1}{\sqrt2}\leq  |x| \leq 1 \\
2\sqrt{1-x^2} - 2\sqrt{\frac{1}{2} - x^2}& \text{ when } 0 \leq  |x| < \frac{1}{\sqrt2} \end{cases} 
$$
does not work, as $l$ isn't Lipschitz at $x=1/\sqrt{2}$, in fact only in $W^{1,s}$ for $s<2$, thus 
$$
Pf(x,y) = \int_{v : (x,v l(x)) \in \Omega} f\left(x,v l(x)\right) dv,
$$
which is indeed a function of $x$, and corresponds to averaging $f$ in the direction $(0,1)$ orthogonal to $\phi$, isn't regular. 
Any hint on how to go for general directions $\phi$ would be great. 
 A: Here is a comment that shows possible obstructions to the existence of such an operator. I'm assuming a bit more regularity here.
Let $(M,g)$ be a compact Riemannian manifold with boundary, and let $\text{vol}_g$ be its Riemannian volume form.
Equip the vector space $C^{\infty}(M)/\mathbb{R}$ of smooth functions integrating to zero with the inner product
$$\langle \alpha,\beta\rangle=\int_{M}g(\text{grad}\,\alpha,\text{grad}\,\beta)\text{vol}_g.$$
For a fixed vector field $\phi$ on M, define the subspace 
$$V_{\phi}=\{f\in C^{\infty}(M)/\mathbb{R}\;\vert\;\text{grad}\,f=\lambda\cdot\phi\text{ for some }\lambda\in C^{\infty}\}.$$
Then I think your problem (correct me if I'm misinterpreting your question) can be reduced to computing for a given $\phi\in\text{Vect}(M)$ and $f\in C^{\infty}(M)/\mathbb{R}$ its orthogonal projection onto $V_{\phi}$ using the inner product above.
Suppose that the projection $Pf$ of $f$ in the above sense exists.
We have
$$0=\int_{M}\lambda g(\text{grad}\,Pf-\text{grad}\,f,\phi)\text{vol}_g,\quad \forall \lambda\in C^{\infty}(M),$$
which implies that $\text{grad}\,Pf-\text{grad}\,f$ has to be orthogonal to $\phi$ at every point in $M$.
Thus, in order to find $Pf$, we can perform orthogonal projection of $\text{grad}\,f$ onto the subspace generated by $\phi$ in each tangent space $T_p M$ using the inner product $g_p$, which leads to a vector field $X_{\phi}$.
We then have to find a potential, i.e. a function $Pf\in C^{\infty}(M)/\mathbb{R}$ such that $\text{grad}\,Pf=X_{\phi}$.
However, the existence of such a function leads to a contradiction if $X_{\phi}=\lambda\cdot\phi$ with some $\lambda\in C^{\infty}$, for which
$$d\lambda\wedge\phi^{\flat}+\lambda d\phi^{\flat}\neq0,$$
where $\phi^{\flat}$ is the one-form dual to $\phi$.

For a counterexample when $M$ is the annulus in $\mathbb{R}^2$ and $g$ is the Euclidean metric on $\mathbb{R}^2$ restricted onto $M$, choose $f=xy$ and $\phi=\partial_x$ (i.e. $\phi=(1,0)$). Then $df=y dx+x dy$ and thus the orthogonal projection onto $\phi$ is given by $X_{\phi}^{\flat}=y dx$. But we have $dX_{\phi}^{\flat}=-\text{vol}_g\neq 0$, so there cannot be a function $Pf$ such that $d(Pf)=X_{\phi}^{\flat}$.
