docs: def centroid_1dg(data, error=None, mask=None):
      """
      Calculate the centroid of a 2D array by fitting 1D Gaussians to the
      marginal ``x`` and ``y`` distributions of the array.
  
      Invalid values (e.g. NaNs or infs) in the ``data`` or ``error``
      arrays are automatically masked.  The mask for invalid values
      represents the combination of the invalid-value masks for the
      ``data`` and ``error`` arrays.
  
      Parameters
      ----------
      data : array_like
          The 2D data array.
  
      error : array_like, optional
          The 2D array of the 1-sigma errors of the input ``data``.
  
      mask : array_like (bool), optional
          A boolean mask, with the same shape as ``data``, where a `True`
          value indicates the corresponding element of ``data`` is masked.
  
      Returns
      -------
      centroid : `~numpy.ndarray`
          The ``x, y`` coordinates of the centroid.
      """
  
      data = np.ma.asanyarray(data)
  
      if mask is not None and mask is not np.ma.nomask:
          mask = np.asanyarray(mask)
          if data.shape != mask.shape:
              raise ValueError('data and mask must have the same shape.')
          data.mask |= mask
  
      if np.any(~np.isfinite(data)):
          data = np.ma.masked_invalid(data)
          warnings.warn('Input data contains input values (e.g. NaNs or infs), '
                        'which were automatically masked.', AstropyUserWarning)
  
      if error is not None:
          error = np.ma.masked_invalid(error)
          if data.shape != error.shape:
              raise ValueError('data and error must have the same shape.')
          data.mask |= error.mask
  
          error.mask = data.mask
          xy_error = np.array([np.sqrt(np.ma.sum(error**2, axis=i))
                               for i in [0, 1]])
          xy_weights = [(1.0 / xy_error[i].clip(min=1.e-30)) for i in [0, 1]]
      else:
          xy_weights = [np.ones(data.shape[i]) for i in [1, 0]]
  
      # assign zero weight to masked pixels
      if data.mask is not np.ma.nomask:
          bad_idx = [np.all(data.mask, axis=i) for i in [0, 1]]
          for i in [0, 1]:
              xy_weights[i][bad_idx[i]] = 0.
  
      xy_data = np.array([np.ma.sum(data, axis=i) for i in [0, 1]])
  
      constant_init = np.ma.min(data)
      centroid = []
      for (data_i, weights_i) in zip(xy_data, xy_weights):
          params_init = gaussian1d_moments(data_i)
          g_init = Const1D(constant_init) + Gaussian1D(*params_init)
          fitter = LevMarLSQFitter()
          x = np.arange(data_i.size)
          g_fit = fitter(g_init, x, data_i, weights=weights_i)
          centroid.append(g_fit.mean_1.value)
  
      return np.array(centroid)