I'm relatively new to Python and I've been trying to obtain contour plots that look similar to this one:
Contour Plot with linear topography
I've more or less successfully done so using the code below:
#---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)) ''' Now I should go on to use real topography data, instead of the ideal linear one h(x) = -x/1000 + 45 used in the code.
Contour Plor with real topography
I have my topography data like so:
x h(x) 0 2 3.5 8 . . . . . . How can I incorporate this information into the code?
Thank you!
