ngrid = 500
x = np.array(lats)
y = np.array(lons)
z = np.array(data)
xi = np.linspace(lat_min, lat_max, ngrid)
yi = np.linspace(lon_min, lon_max, ngrid)
xi, yi = np.meshgrid(xi, yi)
zi = griddata(x, y, z, xi, yi)
# draw the map
plt.xlim(lat_min, lat_max)
plt.ylim(lon_min, lon_max)
plt.contour(xi, yi, zi, 20, linewidths=1)
plt.scatter(x, y, c=z, s=20)