Я относительно новичок в Python, и я пытался получить контурные графики, похожие на этот:
Контурный график с линейной топографией
Я более или менее успешно сделал это, используя приведенный ниже код:
#---Packages---#
import numpy as np
import matplotlib.pyplot as plt
#---Constants and functions---#
import gc
gc.collect()
# Ekman layer depth
D = 60
# Density
rho_0 = 1025
# Coriolis parameter
f = 0.545/(10**4)
# c parameter
c = (1+1j)*np.pi/D
# Zonal geostrophic wind speed
ug = 0
# Wind shear stress
tau_x = 0
tau_y = -0.07
tau = tau_x + 1j*tau_y
# Surface Ekman velocity in the deep ocean limit
u0 = (1-1j)*np.pi*tau/(rho_0*f*D)
# Bathymetry profile
h = lambda x : -x/1000 + 45
zeta = lambda x: np.pi*h(x)/D
# Horizontally-varying structure functions
alpha = lambda x: (np.cosh(zeta(x))*np.cos(zeta(x)))**2 + (np.sinh(zeta(x))*np.sin(zeta(x)))**2
S1 = lambda x: np.cosh(zeta(x))*np.cos(zeta(x))/alpha(x)
S2 = lambda x: np.sinh(zeta(x))*np.sin(zeta(x))/alpha(x)
T1 = lambda x: np.cosh(zeta(x))*np.sinh(zeta(x))/alpha(x)
T2 = lambda x: np.cos(zeta(x))*np.sin(zeta(x))/alpha(x)
# Meridional geostrophic wind speed
vg = lambda x: (2*np.pi*tau_y/(rho_0*f*D))*(1-S1(x))/(T1(x)-T2(x))
# Complex total geostrophic velocity
ugc = lambda x: ug + 1j*vg(x)
# Ekman transport
Uek = np.absolute(tau_y/(rho_0*f))
# Stream function
phi = lambda x, z: np.real((1/c)*(u0*(1-(np.cosh(c*(z+h(x))))/(np.cosh(c*h(x)))) + ugc(x)*np.sinh(c*z)/np.cosh(c*h(x))))
#---Plotting---#
# Meridional geostrophic velocity
'''
of = np.linspace(-100000,0,500)
V_G = vg(of)
plt.plot(of,V_G,'b')
plt.show()
'''
# Contour Plot
x1 = np.linspace(-100000,0,1000)
z1 = np.linspace(-100,0,1000)
breaks = np.linspace(-1, 1, 20)
X1, Z1 = np.meshgrid(x1,z1)
Z = phi(X1,Z1)/Uek
plt.figure()
CS = plt.contour(x1,z1,Z,breaks)
plt.clabel(CS, inline=1, fontsize=10)
plt.grid()
plt.show()
# Value checking
'''
print(h(-40000))
print(zeta(-40000))
print(vg(-40000))
print(Uek)
print(phi(-40000,-300))
'''
Теперь я должен продолжить использовать реальные данные топографии вместо идеального линейного h (x) = -x / 1000 +45 используется в коде.
Contour Plor с реальной топографией
У меня есть данные о моей топографии, например:
x h(x)
0 2
3.5 8
. .
. .
. .
Как включитьэта информация в коде?
Спасибо!