Plotting FiPy equation on a map

Question

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.

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.

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

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

Plotting FiPy equation on a map

0 Answers