Most efficient way to save 2D array and its axes values

Question

Suppose I have some large arrays as a result of a 2-variable function:

import numpy as np

lenx = 100
leny = 37

x = np.linspace(0, 10, lenx)
y = np.linspace(0, 20, leny)
z = np.zeros((lenx, leny))

fx = x + 1 #  Might need this later

for i in range(lenx):
    for k in range(leny):
        z[i, k] = x + y

I can plot this using plt.imshow (for example) after I am done calculating z, but sometimes I want to save this data (x, y, z, and maybe fx too) and replot it later. What would be the best way to go about it?

I could expand this into a 3 by lenx*leny array:

[[x0, y0, z00],
 [x0, y1, z01],
 ...
 [xn, ym, znm]]

But that seems very inefficient (x and y values would be repeated many times unnecessarily) and would take a long time to load every time. I also thought of a json file but my understanding is that they are text-based, so in case of very large arrays (think 200 x 250000) it would take up a lot of memory.

What is the best option in this case?