I would like to solve the following problem $$ \vec{v}(x,y)= k\, \nabla \theta(x,y) $$ with respect to the unknown function $\theta$. Parameter $k$ is just a real constant quantity.
I have two matrices where I have stored the values of $v_x(x,y)$ and $v_y(x,y)$, which are the components of the 2D vector field $\vec{v}$.
You can convert this equation into a Poisson equation and use standard numerical methods to solve it like finite difference to obtain $\theta(x,y)$. If you take a divergence from both side:
$$\nabla^{2} \theta = \frac{1}{k}\nabla \cdot \vec{v}$$
Since, you know the $\vec{v}$ you can approximate its divergence as:
$$\nabla \cdot \vec{v} = \frac{\partial v_{x}}{\partial x} + \frac{\partial v_{y}}{\partial y} = \frac{v_{x}(x+\Delta x,y) - v_{x}(x-\Delta x,y)}{2 \Delta x} + \frac{v_{y}(x,y+\Delta y) - v_{y}(x,y-\Delta y)}{2 \Delta y}$$
The only concern here is that how you want to handle the boundaries. If it makes sense for simplest case you can assume a periodic boundary condition here. Now you have $\nabla \cdot \vec{v}$ stored in an array and for the rest I assume it is known already. In order to find $\theta$ from the defined Poisson equation, you need to know the boundary conditions of $\theta$ but again I assume it's defined as periodic boundary conditions on all boundaries. So, now you can discretize the Poisson equation as:
$$\frac{\theta(x+\Delta x,y)+\theta(x-\Delta x,y)-2\theta(x,y)}{\Delta x^{2}} + \frac{\theta(x,y+\Delta y)+\theta(x,y-\Delta y)-2\theta(x,y)}{\Delta y^{2}} = \frac{1}{k} \mathrm{div}(x,y)$$
Where $\mathrm{div}(x,y)$ is the array that you stores the values of $\nabla \cdot \vec{v}$ in it beforehand. Now you have a standard discretized Poisson equation that you can convert it to a system of linear equations. In fact, I want to convert above system of discretized equations into a matrix form as:
$$A x = b$$
Matrix $A$ is called mass matrix and if you have a 2D grid of $N \times M$ points, $A$'s dimension is $NM \times NM$. Now, let's say I can map $x$ and $y$ grid into matrix coordinate of $i$ and $j$. In fact, you have $NM$ points in your $x$ and $y$ grid and you can assign an ID from 0 to $NM - 1$ to each of them. This would be a function as:
$$F(x,y) = i$$
In fact you want to flatten a 2D grid into 1D as:
$$F(x,y) = x + N y$$
You give $x$ and $y$ coordinate and it just gives you the id of the point. Then I can assemble $A$ and $b$ as:
$$A_{F(x,y),F(x,y)} = -2 \Bigg ( \frac{1}{\Delta x^{2}} + \frac{1}{\Delta y^{2}} \Bigg)$$
$$A_{F(x+\Delta x,y),F(x,y)} = \frac{1}{\Delta x^{2}}$$
$$A_{F(x-\Delta x,y),F(x,y)} = \frac{1}{\Delta x^{2}}$$
$$A_{F(x,y + \Delta y),F(x,y)} = \frac{1}{\Delta y^{2}}$$
$$A_{F(x,y - \Delta y),F(x,y)} = \frac{1}{\Delta y^{2}}$$
$$b_{F(x,y)} = \frac{1}{k} \mathrm{div}(x,y)$$
Keep in mind you need to treat boundaries specially based on your boundary condition to assemble $A$ and $b$. Finally you solve this matrix equation with standard matrix solvers available in Python or any other programming language that you use for your problem.
