ValueError plotting contourf in Cartopy Azimuthal Equidistant map projection

I would like to use contourf in a Cartopy Azimuth Equidistant map projection. For context, I am trying to plot travel time (in h) of a signal across the globe. Roughly, this is somewhat what I am trying my plot to look like:

(Image credits to @htonchia, https://github.com/SciTools/cartopy/issues/1421)

However, when I try to plot it, it gives the error: line 185, in geos_multipolygon_from_polygons ValueError: Sequences of multi-polygons are not valid arguments

To reproduce:

# Data
# Longitudes of stations
longs = [-171.7827,  -171.7827,   179.1966,   179.1966,  -159.7733,  -159.7733,
  174.7043,   174.7043,   172.9229,   172.9229,   159.9475,   159.9475,
 -157.4457, -157.4457,   146.24998,  146.24998, -169.5292,  -169.5292,
  166.652,    166.652,   -155.5326,  -155.5326,  -158.0112,  -158.0112,
 -177.3698,  -177.3698,   144.8684,   166.7572,   166.7572,   117.239,
  117.239,    125.5791,   125.5791,   110.5354,   110.5354,   144.4382,
  144.4382,   138.20406,  138.20406, -176.6842,  -176.6842,   121.4971,
  121.4971,   126.62436,  126.62436,  -64.0489,   -64.0489,  -123.3046,
 -123.3046,  -110.7847,  -110.7847,   -90.2861,   -90.2861,  -106.4572,
 -106.4572,  -106.4572,  -147.8616,  -147.8616,  -147.8616,  -104.0359,
 -104.0359,   -95.83812,  -95.83812,  -70.7005,   -70.7005,    98.9443,
   98.9443,   -88.2763,   -88.2763,   -61.9787,  -61.9787,   -78.4508,
  -78.4508,  -175.385  ]

# Latitudes of stations
lats = [-13.9085,   -13.9085,    -8.5259,    -8.5259,   -21.2125,   -21.2125,
 -41.3087,   -41.3087,     1.3549,     1.3549,    -9.4387,    -9.4387,
   2.0448,     2.0448,   -20.08765,  -20.08765,   16.7329,    16.7329,
  19.2834,    19.2834,    19.7573,    19.7573,    21.42,      21.42,
  28.2156,    28.2156,    13.5893,   -77.8492,   -77.8492,   -32.9277,
 -32.9277,     7.0697,     7.0697,   -66.2792,   -66.2792,   -89.9289,
 -89.9289,    36.54567,   36.54567,   51.8823,    51.8823,    24.9735,
  24.9735,    37.47768,   37.47768, -64.7744,   -64.7744,    44.5855,
  44.5855,    32.3098,    32.3098,    -0.6742,    -0.6742,    34.94591,
  34.94591,   34.94591,   64.873599,  64.873599,  64.873599,  44.1212,
  44.1212,    29.96478,   29.96478,  -29.011,    -29.011,     18.8141,
  18.8141,    20.2263,    20.2263,   -38.0568,   -38.0568,     0.2376,
   0.2376,   -20.57    ]

# Time (h) signal detected after eruption
travel_time_h = [ 0.95296297,  0.95332528,  1.49046297,  1.4905475,   1.67046297, 
1.67026972, 2.3705475,   2.37046297,  2.60249194,  2.60240741,  2.7537963,   2.75360306,
  3.00943639,  3.00935186,  3.65610306,  3.65601852,  3.93165861,  3.93157408,
 16.13526972,  4.43074074,  4.61268519,  4.6130475,   4.6730475,   4.67296297,
  5.01026972,  5.01046297,  5.20768519,  5.96546297,  5.9655475,   6.49693639,
  6.49685186,  6.40324074,  6.40332528,  6.53740741,  6.53721417,  7.12074074,
  7.1205475,   7.34546297,  7.34499194,  7.26157408,  7.26221417,  7.64546297,
  7.6455475,   8.13407408,  8.13388083,  7.97693639,  7.97740741,  8.05082528,
  8.05101852,  8.00240741,  8.00221417,  8.65943639,  8.65907408,  8.41907408,
  8.41776972,  8.42722222,  8.94324074,  8.9430475,   8.94333333,  9.2555475,
  9.25601852,  8.99240741,  8.99249194,  9.26851852,  9.26749194,  9.16165861,
  9.16185186,  9.41990741,  9.41999194,  9.30851852,  9.31360306,  9.82324074,
  9.82332528,  0.        ]

I then interpolate the data to use in contourf

import matplotlib.pyplot as plt
import numpy as np
from cartopy import crs as ccrs
from scipy.interpolate import griddata

# Interpolate for contour
X, Y = np.meshgrid(longs, lats)
Z = griddata((longs, lats), travel_time_h, (X, Y), method='linear')

And try to plot it using Cartopy:

# Initialize figure
fig = plt.figure(figsize=(10, 8))
projLae = ccrs.AzimuthalEquidistant(central_longitude=-175.385, central_latitude=-20.57)
ax = plt.subplot(1, 1, 1, projection=projLae)

# Plot contour first as background
start_h, end_h, interval_h = 0.0, 10.0, 0.5
levels = np.arange(start=start_h, stop=end_h, step=interval_h)  # levels of contour
contour = ax.contourf(X, Y, Z, levels=levels, vmin=start_h, vmax=end_h, transform=ccrs.PlateCarree())

# Add colorbar for contour
cbar = fig.colorbar(contour, orientation='horizontal')   
cbar.ax.set_xlabel(f"Time [hr]")

# Plot station locations
ax.scatter(longs, lats, s=8, marker='*', color='red', transform=ccrs.PlateCarree())

# Plot map details
ax.coastlines()
ax.set_global()
plt.show()

Not sure what is going on, if is a ax.contourf problem and/or a Cartopy Azimuthal Equidistant projection problem. I'm using Cartopy version 0.21.1.

I appreciate any help!

There are 2 major steps that need amendment. Firstly, the creation of gridmesh for use in data interpolation, and the data interpolation step. See comments in the code for details.

Updated code:

import matplotlib.pyplot as plt
import numpy as np
from cartopy import crs as ccrs
from scipy.interpolate import griddata

# Data
# Longitudes of stations
longs = [-171.7827,  -171.7827,   179.1966,   179.1966,  -159.7733,  -159.7733,
  174.7043,   174.7043,   172.9229,   172.9229,   159.9475,   159.9475,
 -157.4457, -157.4457,   146.24998,  146.24998, -169.5292,  -169.5292,
  166.652,    166.652,   -155.5326,  -155.5326,  -158.0112,  -158.0112,
 -177.3698,  -177.3698,   144.8684,   166.7572,   166.7572,   117.239,
  117.239,    125.5791,   125.5791,   110.5354,   110.5354,   144.4382,
  144.4382,   138.20406,  138.20406, -176.6842,  -176.6842,   121.4971,
  121.4971,   126.62436,  126.62436,  -64.0489,   -64.0489,  -123.3046,
 -123.3046,  -110.7847,  -110.7847,   -90.2861,   -90.2861,  -106.4572,
 -106.4572,  -106.4572,  -147.8616,  -147.8616,  -147.8616,  -104.0359,
 -104.0359,   -95.83812,  -95.83812,  -70.7005,   -70.7005,    98.9443,
   98.9443,   -88.2763,   -88.2763,   -61.9787,  -61.9787,   -78.4508,
  -78.4508,  -175.385  ]

# Latitudes of stations
lats = [-13.9085,   -13.9085,    -8.5259,    -8.5259,   -21.2125,   -21.2125,
 -41.3087,   -41.3087,     1.3549,     1.3549,    -9.4387,    -9.4387,
   2.0448,     2.0448,   -20.08765,  -20.08765,   16.7329,    16.7329,
  19.2834,    19.2834,    19.7573,    19.7573,    21.42,      21.42,
  28.2156,    28.2156,    13.5893,   -77.8492,   -77.8492,   -32.9277,
 -32.9277,     7.0697,     7.0697,   -66.2792,   -66.2792,   -89.9289,
 -89.9289,    36.54567,   36.54567,   51.8823,    51.8823,    24.9735,
  24.9735,    37.47768,   37.47768, -64.7744,   -64.7744,    44.5855,
  44.5855,    32.3098,    32.3098,    -0.6742,    -0.6742,    34.94591,
  34.94591,   34.94591,   64.873599,  64.873599,  64.873599,  44.1212,
  44.1212,    29.96478,   29.96478,  -29.011,    -29.011,     18.8141,
  18.8141,    20.2263,    20.2263,   -38.0568,   -38.0568,     0.2376,
   0.2376,   -20.57    ]

# Time (h) signal detected after eruption
travel_time_h = [ 0.95296297,  0.95332528,  1.49046297,  1.4905475,   1.67046297, 
1.67026972, 2.3705475,   2.37046297,  2.60249194,  2.60240741,  2.7537963,   2.75360306,
  3.00943639,  3.00935186,  3.65610306,  3.65601852,  3.93165861,  3.93157408,
 16.13526972,  4.43074074,  4.61268519,  4.6130475,   4.6730475,   4.67296297,
  5.01026972,  5.01046297,  5.20768519,  5.96546297,  5.9655475,   6.49693639,
  6.49685186,  6.40324074,  6.40332528,  6.53740741,  6.53721417,  7.12074074,
  7.1205475,   7.34546297,  7.34499194,  7.26157408,  7.26221417,  7.64546297,
  7.6455475,   8.13407408,  8.13388083,  7.97693639,  7.97740741,  8.05082528,
  8.05101852,  8.00240741,  8.00221417,  8.65943639,  8.65907408,  8.41907408,
  8.41776972,  8.42722222,  8.94324074,  8.9430475,   8.94333333,  9.2555475,
  9.25601852,  8.99240741,  8.99249194,  9.26851852,  9.26749194,  9.16165861,
  9.16185186,  9.41990741,  9.41999194,  9.30851852,  9.31360306,  9.82324074,
  9.82332528,  0.        ]

# Create ranged values for making gridmesh
lon_min, lon_max, lat_min, lat_max = min(longs), max(longs), min(lats), max(lats)
grid_x, grid_y = 32, 24   # grid columns, rows
alons = np.linspace(lon_min, lon_max, grid_x)
alats = np.linspace(lat_min, lat_max, grid_y)

# Original bad code: X, Y = np.meshgrid(longs, lats)
X, Y = np.meshgrid(alons , alats)

# Interpolating for Z, need proper `fill_value`.
# - Available data points: longs, lats, travel_time_h
# - Places to interpolate for data: (X, Y)
# If interpolation fails, `fill_value` is used.
Z = griddata((longs, lats), travel_time_h, (X, Y), method='linear', fill_value=0.01) #linear, cubic

# Set projection
projLae = ccrs.AzimuthalEquidistant(central_longitude=-175.385, central_latitude=-20.57)

# Set figure, axes for plotting
fig = plt.figure(figsize = (10,8))
ax = fig.add_subplot(111, projection=projLae)

# Plot contour
start_h, end_h, interval_h = 0.0, 10.0, 1
levels = np.arange(start=start_h, stop=end_h, step=interval_h)  # levels of contour

contour = ax.contourf(X, Y, Z, levels=levels, vmin=start_h, vmax=end_h, transform=ccrs.PlateCarree(), zorder=20, alpha=0.75)

# Add colorbar for contour
cbar = fig.colorbar(contour, orientation='horizontal', shrink=0.45)   
cbar.ax.set_xlabel(f"Time [hr]")

# Plot station locations
ax.scatter(longs, lats, s=8, marker='*', color='red', transform=ccrs.PlateCarree(), zorder=30)

# Plot map details
ax.coastlines()
ax.set_global()

plt.show()

Output plot:

Unfortunately, the plot has small white gap at the longitude of the dateline. This white gap can be eliminated by adding cyclic data and regenerate the contour. The relevant code and the final plot are as follows.

# Add cyclic data to eliminate the white gap
import cartopy.util as cutil
cdata, clon, clat = cutil.add_cyclic(Z, X, Y)
contour = ax.contourf(clon, clat, cdata, levels=levels, vmin=start_h, vmax=end_h, transform=ccrs.PlateCarree(), zorder=20, alpha=0.75)