docs: def fit_2dgaussian(data, error=None, mask=None):
      """
      Fit a 2D Gaussian plus a constant to a 2D image.
  
      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 array of the image.
  
      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
      -------
      result : A `GaussianConst2D` model instance.
          The best-fitting Gaussian 2D model.
      """
  
      from ..morphology import data_properties  # prevent circular imports
  
      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
          weights = 1.0 / error.clip(min=1.e-30)
      else:
          weights = np.ones(data.shape)
  
      if np.ma.count(data) < 7:
          raise ValueError('Input data must have a least 7 unmasked values to '
                           'fit a 2D Gaussian plus a constant.')
  
      # assign zero weight to masked pixels
      if data.mask is not np.ma.nomask:
          weights[data.mask] = 0.
  
      mask = data.mask
      data.fill_value = 0.0
      data = data.filled()
  
      # Subtract the minimum of the data as a crude background estimate.
      # This will also make the data values positive, preventing issues with
      # the moment estimation in data_properties (moments from negative data
      # values can yield undefined Gaussian parameters, e.g. x/y_stddev).
      props = data_properties(data - np.min(data), mask=mask)
  
      init_const = 0.    # subtracted data minimum above
      init_amplitude = np.ptp(data)
      g_init = GaussianConst2D(constant=init_const, amplitude=init_amplitude,
                               x_mean=props.xcentroid.value,
                               y_mean=props.ycentroid.value,
                               x_stddev=props.semimajor_axis_sigma.value,
                               y_stddev=props.semiminor_axis_sigma.value,
                               theta=props.orientation.value)
      fitter = LevMarLSQFitter()
      y, x = np.indices(data.shape)
      gfit = fitter(g_init, x, y, data, weights=weights)
  
      return gfit