landarea

# my answer
from tqdm import tqdm

subaerial = np.zeros(golfcube.shape[0])
topset = np.zeros(golfcube.shape[0])
for T in tqdm(range(golfcube.shape[0])):

  # do it with just pixels above sea level
  Land_A = golfcube['eta'][T,:,:] > golfcube.meta['H_SL'][T]
  Land_C = np.count_nonzero(Land_A)
  subaerial[T] = Land_C * golfcube.meta["dx"] * golfcube.meta["dx"]

  # do it with Opening Angle method
  lm = spl.mask.LandMask(
      golfcube['eta'][T, :, :],
      elevation_threshold=golfcube.meta['H_SL'][T])
  topset[T] = np.sum(lm.mask) * golfcube.meta["dx"] * golfcube.meta["dx"]

fig, ax = plt.subplots()

ax.plot(golfcube.t, subaerial / 1e6, label="subaerial")
ax.plot(golfcube.t, topset / 1e6, label="topset")

ax.set_ylabel("area [km$^2$]")
ax.set_xlabel("time [elapsed seconds]")
ax.legend()
plt.show()