0

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!

ASPVL
  • 1
  • 1
  • As I understand it, the question is more about how to import data to a numpy array than plotting a contour graph?, see for example https://stackoverflow.com/q/4315506/8069403 – xdze2 Oct 23 '19 at 13:41
  • @xdze2 Thank you for your answer. Actually I'm interested in how to make the contour plot from that array. Everything in the code depends on the function h(x), but in this case, there's not a function, only values of that function for particular values of x. – ASPVL Oct 24 '19 at 10:53
  • I think the function to use is [tricontour](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.tricontour.html), which works on unstructured data (contrary to contour which need data on a grid). I would be helpfull if you rewrite the question toward this points, and simplify the example code including example data which could be copy/pasted (see [mcve]) – xdze2 Oct 24 '19 at 11:18

0 Answers0