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.

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$$

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.