5

I would like to model the biodiffusion of stork through migration.

I have an equation that gives me that kind of results, and I would like to plot it on a map to represent the movements.

enter image description here

My code equation :

D = Variable(value=(0))
ux = Variable(value=0)
uy = Variable(value=0)
q  = 0


Smax = 1.5


Lx, Ly = 2500, 2500
Nx, Ny = 100, 100
dx, dy = Lx/Nx, Ly/Ny
Cx, Cy = 1.25*Lx/4, 3*Ly/4
sigma  = Lx/10
T  = 10
dt = T/1000
Np = 10     # Number of snapshots of the solution (one every Nt/Np time steps)
Nt = Np*np.ceil(T/(Np*dt)).astype('int')   # Nt must be an integer divisible by Np



# Define the grid/mesh
mesh = Grid2D(nx=Nx, ny=Ny, dx=dx, dy=dy)
x, y = mesh.cellCenters[0], mesh.cellCenters[1]

# Define the model variable and set the boundary conditions
phi = CellVariable(name="numerical solution", mesh=mesh, value=Smax*numerix.exp(-((x-Cx)**2 + (y-Cy)**2)/(0.5*sigma**2)))
meshBnd = mesh.facesLeft | mesh.facesRight | mesh.facesBottom | mesh.facesTop




# impose zero flux on all 4 boundaries
phi.faceGrad.constrain(0., where=meshBnd)

# Define and then solve the equation
eq = TransientTerm() == DiffusionTerm(coeff=D) \
    - ExponentialConvectionTerm(coeff=(ux,uy)) \
    + ImplicitSourceTerm(coeff=q)
my_sol = np.zeros((Np,Nx*Ny))      # Matrix with Np solution snapshots
my_sol[0,:] = phi
k = 1

for step in np.arange(1,Nt):
    if step < Nt/2:
        print(step)
        D.setValue(0)

    else:
        D.setValue(5000)
        ux.setValue(100)
        uy.setValue(-150)
    eq.solve(var=phi, dt=dt)
    if np.mod(step,Nt/Np)==0:
        print(step,k)
        my_sol[k,:] = phi
        k += 1


# Plot the solution
fig, ax = plt.subplots(1,1)
img = plt.imread("map.png")

ax.imshow(img, extent=[0, 2500, 0, 2500])


#Viewer(vars=phi, datamin=0., datamax=100.)
xg, yg = np.meshgrid(np.linspace(0,Lx,Nx+1), np.linspace(0,Ly,Ny+1))
xd, yd = np.meshgrid(np.linspace(dx/2,Lx,Nx), np.linspace(dy/2,Ly,Ny))
for i in np.arange(Np):
    sol = my_sol[i,:].reshape((Nx,Ny))
    sol = griddata((xd.ravel(),yd.ravel()), my_sol[i,:], (xg, yg), method='nearest')
    plt.contourf(xg, yg, sol, zorder=0)
    plt.clim(0,1)
    ax.set_aspect('equal', 'box')
    #plt.savefig('./Plots/fipy_advdiff_'+str(i)+'.png', bbox_inches='tight',pad_inches=0.0)
    plt.imshow(img, zorder=1, extent=[0, 2500, 0, 2500])
    plt.show()

I get this as result when I put zorder=1 for the image so it works, but when I put zorder=0 I can't see the map because the background of my solution (shown above) is not transparent.

enter image description here

@tandat suggest to had an alpha coefficient to my image to make it more transparent : this is the result with alpha = 0.2 :

enter image description here

I am already satisfied with the result but I don't close this question yet in case anyone gets another idea !

Jouline
  • 90
  • 7
  • 1
    As I can not run the code so I'll make a quick comment instead of a full answer. You can try setting some alpha to make images overlay. Syntax will look like `plt.imshow(..., alpha=...)` – tandat Dec 15 '21 at 13:18
  • @tandat What do I need to add for you to be able to run it ? Here are the requirements : `from fipy import CellVariable, FaceVariable, Grid2D, ExponentialConvectionTerm, TransientTerm, DiffusionTerm, ImplicitSourceTerm, boundaryConditions, FixedFlux, Viewer from fipy.tools import numerix from fipy.variables.variable import Variable import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import griddata` Btw it's working with your idea ! I'll edit my question! – Jouline Dec 15 '21 at 13:42
  • 2
    Try setting your background image with `zorder=0, alpha=1` because it is background. And overlay the function graph on the image with `zorder=1, alpha=0.2` (or some smaller `alpha`). It can help the plot look better, I guess. And if my comment solved your answer, please mark it as `solved`. – tandat Dec 15 '21 at 13:50
  • 2
    @tandat I edited my question if you want to look at the result, I used the exact same coefficient for `alpha`. I just want to know if there is another answer to this question, but I will close it if no one else response within 24 hours. – Jouline Dec 15 '21 at 13:56
  • I'd be inclined to try a [sequential colormap](https://matplotlib.org/stable/tutorials/colors/colormaps.html#sequential) and overplot the data on the map, e.g., plot phi with `zorder=1, alpha=.2, cmap="Purples"` and the map with `zorder=0, alpha=1`. If @tandat wrote up a formal answer, I would upvote it. – jeguyer Dec 15 '21 at 15:11

0 Answers0