Just use numpy logical and there where statement to apply a 1D interpolation.

import numpy as np
from scipy import interpolate

def fill_nan(A):
    '''
    interpolate to fill nan values
    '''
    inds = np.arange(A.shape[0])
    good = np.where(np.isfinite(A))
    f = interpolate.interp1d(inds[good], A[good],bounds_error=False)
    B = np.where(np.isfinite(A),A,f(inds))
    return B

Slightly optimized version based on response of BRYAN WOODS. It handles starting and ending values of source data correctly, and it is faster on 25-30% than original version. Also you may use different kinds of interpolations (see scipy.interpolate.interp1d documentations for details).

import numpy as np
from scipy.interpolate import interp1d

def fill_nans_scipy1(padata, pkind='linear'):
"""
Interpolates data to fill nan values

Parameters:
    padata : nd array 
        source data with np.NaN values
    
Returns:
    nd array 
        resulting data with interpolated values instead of nans
"""
aindexes = np.arange(padata.shape[0])
agood_indexes, = np.where(np.isfinite(padata))
f = interp1d(agood_indexes
           , padata[agood_indexes]
           , bounds_error=False
           , copy=False
           , fill_value="extrapolate"
           , kind=pkind)
return f(aindexes)

In [17]: adata = np.array([1, 2, np.NaN, 4])
Out[18]: array([ 1.,  2., nan,  4.])
In [19]: fill_nans_scipy1(adata)
Out[19]: array([1., 2., 3., 4.])

I use the interpolation for replacing all NaN values.

A = np.array([1, nan, nan, 2, 2, nan, 0])
np.interp(np.arange(len(A)), 
          np.arange(len(A))[np.isnan(A) == False], 
          A[np.isnan(A) == False])

Output :

array([1. , 1.33333333, 1.66666667, 2. , 2. , 1. , 0. ])

Interpolation and extrapolation with padding keywords

The following solution interpolates the nan values in an array by np.interp, if a finite value is present on both sides. Nan values at the borders are handled by np.pad with modes like constant or reflect.

    import numpy as np
    import matplotlib.pyplot as plt
    
    
    def extrainterpolate_nans_1d(
            arr, kws_pad=({'mode': 'edge'}, {'mode': 'edge'})
            ):
        """Interpolates and extrapolates nan values.
    
        Interpolation is linear, compare np.interp(..).
        Extrapolation works with pad keywords, compare np.pad(..).
    
        Parameters
        ----------
        arr : np.ndarray, shape (N,)
            Array to replace nans in.
        kws_pad : dict or (dict, dict)
            kwargs for np.pad on left and right side
    
        Returns
        -------
        bool
            Description of return value
    
        See Also
        --------
        https://numpy.org/doc/stable/reference/generated/numpy.interp.html
        https://numpy.org/doc/stable/reference/generated/numpy.pad.html
        https://stackoverflow.com/a/43821453/7128154
        """
        assert arr.ndim == 1
        if isinstance(kws_pad, dict):
            kws_pad_left = kws_pad
            kws_pad_right = kws_pad
        else:
            assert len(kws_pad) == 2
            assert isinstance(kws_pad[0], dict)
            assert isinstance(kws_pad[1], dict)
            kws_pad_left = kws_pad[0]
            kws_pad_right = kws_pad[1]
    
        arr_ip = arr.copy()
    
        # interpolation
        inds = np.arange(len(arr_ip))
        nan_msk = np.isnan(arr_ip)
        arr_ip[nan_msk] = np.interp(inds[nan_msk], inds[~nan_msk], arr[~nan_msk])
    
        # detemine pad range
        i0 = next(
            (ids for ids, val in np.ndenumerate(arr) if not np.isnan(val)), 0)[0]
        i1 = next(
            (ids for ids, val in np.ndenumerate(arr[::-1]) if not np.isnan(val)), 0)[0]
        i1 = len(arr) - i1
        # print('pad in range [0:{:}] and [{:}:{:}]'.format(i0, i1, len(arr)))
    
        # pad
        arr_pad = np.pad(
            arr_ip[i0:], pad_width=[(i0, 0)], **kws_pad_left)
        arr_pad = np.pad(
            arr_pad[:i1], pad_width=[(0, len(arr) - i1)], **kws_pad_right)
    
        return arr_pad
    
    
    # setup data
    ys = np.arange(30, dtype=float)**2/20
    ys[:5] = np.nan
    ys[20:] = 20
    ys[28:] = np.nan
    ys[[7, 13, 14, 18, 22]] = np.nan
    
    
    ys_ie0 = extrainterpolate_nans_1d(ys)
    kws_pad_sym = {'mode': 'symmetric'}
    kws_pad_const7 = {'mode': 'constant', 'constant_values':7.}
    ys_ie1 = extrainterpolate_nans_1d(ys, kws_pad=(kws_pad_sym, kws_pad_const7))
    ys_ie2 = extrainterpolate_nans_1d(ys, kws_pad=(kws_pad_const7, kws_pad_sym))
    
    fig, ax = plt.subplots()
    
    
    ax.scatter(np.arange(len(ys)), ys, s=15**2, label='ys')
    ax.scatter(np.arange(len(ys)), ys_ie0, s=8**2, label='ys_ie0, left_pad edge, right_pad edge')
    ax.scatter(np.arange(len(ys)), ys_ie1, s=6**2, label='ys_ie1, left_pad symmetric, right_pad 7')
    ax.scatter(np.arange(len(ys)), ys_ie2, s=4**2, label='ys_ie2, left_pad 7, right_pad symmetric')
    ax.legend()

As suggested by an earlier comment, the best way to do this is to use a peer reviewed implementation. The pandas library has an interpolation method for 1d data, which interpolates np.nan values in Series or DataFrame:

pandas.Series.interpolate or pandas.DataFrame.interpolate

The documentation is very concise, recommend reading through! My implementation:

import pandas as pd

magnitudes_series = pd.Series(magnitudes)    # Convert np.array to pd.Series
magnitudes_series.interpolate(
    # I used "akima" because the second derivative of my data has frequent drops to 0
    method=interpolation_method,

    # Interpolate from both sides of the sequence, up to you (made sense for my data)
    limit_direction="both",

    # Interpolate only np.nan sequences that have number sequences at the ends of the respective np.nan sequences
    limit_area="inside",

    inplace=True,
)

# I chose to remove np.nan at the tails of data sequence
magnitudes_series.dropna(inplace=True)

result_in_numpy_array = magnitudes_series.values

Importing scipy looks like overkill to me. Here's a simple way using numpy and maintaining the same conventions as np.interp

   def interp_nans(x:[float],left=None, right=None, period=None)->[float]:
    """ 
      e.g. [1 1 1 nan nan 2 2 nan 0] -> [1 1 1 1.3 1.6 2 2  1  0]
    
    """
    xp = [i for i, yi in enumerate(x) if np.isfinite(yi)]
    fp = [yi for i, yi in enumerate(x) if np.isfinite(yi)]
    return list(np.interp(x=list(range(len(x))), xp=xp, fp=fp,left=left,right=right,period=period))