
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples/example_two_dimensional_peak.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_example_two_dimensional_peak.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_example_two_dimensional_peak.py:


Fit Two Dimensional Peaks
=========================

This example illustrates how to handle two-dimensional data with lmfit.

.. GENERATED FROM PYTHON SOURCE LINES 8-15

.. code-block:: Python

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

    import lmfit
    from lmfit.lineshapes import gaussian2d, lorentzian








.. GENERATED FROM PYTHON SOURCE LINES 16-21

Two-dimensional Gaussian
------------------------
We start by considering a simple two-dimensional gaussian function, which
depends on coordinates `(x, y)`. The most general case of experimental
data will be irregularly sampled and noisy. Let's simulate some:

.. GENERATED FROM PYTHON SOURCE LINES 21-28

.. code-block:: Python

    npoints = 10000
    np.random.seed(2021)
    x = np.random.rand(npoints)*10 - 4
    y = np.random.rand(npoints)*5 - 3
    z = gaussian2d(x, y, amplitude=30, centerx=2, centery=-.5, sigmax=.6, sigmay=.8)
    z += 2*(np.random.rand(*z.shape)-.5)
    error = np.sqrt(z+1)







.. GENERATED FROM PYTHON SOURCE LINES 29-30

To plot this, we can interpolate the data onto a grid.

.. GENERATED FROM PYTHON SOURCE LINES 30-40

.. code-block:: Python

    X, Y = np.meshgrid(np.linspace(x.min(), x.max(), 100),
                       np.linspace(y.min(), y.max(), 100))
    Z = griddata((x, y), z, (X, Y), method='linear', fill_value=0)

    fig, ax = plt.subplots()
    art = ax.pcolor(X, Y, Z, shading='auto')
    plt.colorbar(art, ax=ax, label='z')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.show()



.. image-sg:: /examples/images/sphx_glr_example_two_dimensional_peak_001.png
   :alt: example two dimensional peak
   :srcset: /examples/images/sphx_glr_example_two_dimensional_peak_001.png, /examples/images/sphx_glr_example_two_dimensional_peak_001_3_00x.png 3.00x
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 41-42

In this case, we can use a built-in model to fit

.. GENERATED FROM PYTHON SOURCE LINES 42-46

.. code-block:: Python

    model = lmfit.models.Gaussian2dModel()
    params = model.guess(z, x, y)
    result = model.fit(z, x=x, y=y, params=params, weights=1/error)
    lmfit.report_fit(result)




.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    [[Fit Statistics]]
        # fitting method   = leastsq
        # function evals   = 87
        # data points      = 10000
        # variables        = 5
        chi-square         = 20618.1774
        reduced chi-square = 2.06284916
        Akaike info crit   = 7245.87992
        Bayesian info crit = 7281.93162
        R-squared          = 0.87508800
    [[Variables]]
        amplitude:  27.4195833 +/- 0.65062974 (2.37%) (init = 16.51399)
        centerx:    1.99705425 +/- 0.01405864 (0.70%) (init = 1.940764)
        centery:   -0.49516158 +/- 0.01907800 (3.85%) (init = -0.5178641)
        sigmax:     0.54740777 +/- 0.01224965 (2.24%) (init = 1.666582)
        sigmay:     0.73300589 +/- 0.01617042 (2.21%) (init = 0.8332836)
        fwhmx:      1.28904676 +/- 0.02884573 (2.24%) == '2.3548200*sigmax'
        fwhmy:      1.72609693 +/- 0.03807842 (2.21%) == '2.3548200*sigmay'
        height:     10.8758308 +/- 0.35409534 (3.26%) == '0.1591549*amplitude/(max(1e-15, sigmax)*max(1e-15, sigmay))'
    [[Correlations]] (unreported correlations are < 0.100)
        C(amplitude, sigmay) = +0.2464
        C(amplitude, sigmax) = +0.2314




.. GENERATED FROM PYTHON SOURCE LINES 47-49

To check the fit, we can evaluate the function on the same grid we used
before and make plots of the data, the fit and the difference between the two.

.. GENERATED FROM PYTHON SOURCE LINES 49-77

.. code-block:: Python

    fig, axs = plt.subplots(2, 2, figsize=(10, 10))

    vmax = np.nanpercentile(Z, 99.9)

    ax = axs[0, 0]
    art = ax.pcolor(X, Y, Z, vmin=0, vmax=vmax, shading='auto')
    plt.colorbar(art, ax=ax, label='z')
    ax.set_title('Data')

    ax = axs[0, 1]
    fit = model.func(X, Y, **result.best_values)
    art = ax.pcolor(X, Y, fit, vmin=0, vmax=vmax, shading='auto')
    plt.colorbar(art, ax=ax, label='z')
    ax.set_title('Fit')

    ax = axs[1, 0]
    fit = model.func(X, Y, **result.best_values)
    art = ax.pcolor(X, Y, Z-fit, vmin=0, vmax=10, shading='auto')
    plt.colorbar(art, ax=ax, label='z')
    ax.set_title('Data - Fit')

    for ax in axs.ravel():
        ax.set_xlabel('x')
        ax.set_ylabel('y')
    axs[1, 1].remove()
    plt.show()





.. image-sg:: /examples/images/sphx_glr_example_two_dimensional_peak_002.png
   :alt: Data, Fit, Data - Fit
   :srcset: /examples/images/sphx_glr_example_two_dimensional_peak_002.png, /examples/images/sphx_glr_example_two_dimensional_peak_002_3_00x.png 3.00x
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 78-84

Two-dimensional off-axis Lorentzian
-----------------------------------
We now go on to show a harder example, in which the peak has a Lorentzian
profile and an off-axis anisotropic shape. This can be handled by applying a
suitable coordinate transform and then using the `lorentzian` function that
lmfit provides in the lineshapes module.

.. GENERATED FROM PYTHON SOURCE LINES 84-100

.. code-block:: Python

    def lorentzian2d(x, y, amplitude=1., centerx=0., centery=0., sigmax=1., sigmay=1.,
                     rotation=0):
        """Return a two dimensional lorentzian.

        The maximum of the peak occurs at ``centerx`` and ``centery``
        with widths ``sigmax`` and ``sigmay`` in the x and y directions
        respectively. The peak can be rotated by choosing the value of ``rotation``
        in radians.
        """
        xp = (x - centerx)*np.cos(rotation) - (y - centery)*np.sin(rotation)
        yp = (x - centerx)*np.sin(rotation) + (y - centery)*np.cos(rotation)
        R = (xp/sigmax)**2 + (yp/sigmay)**2

        return 2*amplitude*lorentzian(R)/(np.pi*sigmax*sigmay)









.. GENERATED FROM PYTHON SOURCE LINES 101-102

Data can be simulated and plotted in the same way as we did before.

.. GENERATED FROM PYTHON SOURCE LINES 102-121

.. code-block:: Python

    npoints = 10000
    x = np.random.rand(npoints)*10 - 4
    y = np.random.rand(npoints)*5 - 3
    z = lorentzian2d(x, y, amplitude=30, centerx=2, centery=-.5, sigmax=.6,
                     sigmay=1.2, rotation=30*np.pi/180)
    z += 2*(np.random.rand(*z.shape)-.5)
    error = np.sqrt(z+1)

    X, Y = np.meshgrid(np.linspace(x.min(), x.max(), 100),
                       np.linspace(y.min(), y.max(), 100))
    Z = griddata((x, y), z, (X, Y), method='linear', fill_value=0)

    fig, ax = plt.subplots()
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    art = ax.pcolor(X, Y, Z, shading='auto')
    plt.colorbar(art, ax=ax, label='z')
    plt.show()




.. image-sg:: /examples/images/sphx_glr_example_two_dimensional_peak_003.png
   :alt: example two dimensional peak
   :srcset: /examples/images/sphx_glr_example_two_dimensional_peak_003.png, /examples/images/sphx_glr_example_two_dimensional_peak_003_3_00x.png 3.00x
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 122-128

To fit, create a model from the function. Don't forget to tell lmfit that both
`x` and `y` are independent variables. Keep in mind that lmfit will take the
function keywords as default initial guesses in this case and that it will not
know that certain parameters only make physical sense over restricted ranges.
For example, peak widths should be positive and the rotation can be restricted
over a quarter circle.

.. GENERATED FROM PYTHON SOURCE LINES 128-138

.. code-block:: Python

    model = lmfit.Model(lorentzian2d, independent_vars=['x', 'y'])
    params = model.make_params(amplitude=10, centerx=x[np.argmax(z)],
                               centery=y[np.argmax(z)])
    params['rotation'].set(value=.1, min=0, max=np.pi/2)
    params['sigmax'].set(value=1, min=0)
    params['sigmay'].set(value=2, min=0)

    result = model.fit(z, x=x, y=y, params=params, weights=1/error)
    lmfit.report_fit(result)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    [[Fit Statistics]]
        # fitting method   = leastsq
        # function evals   = 73
        # data points      = 10000
        # variables        = 6
        chi-square         = 11287.3823
        reduced chi-square = 1.12941588
        Akaike info crit   = 1223.00402
        Bayesian info crit = 1266.26606
        R-squared          = 0.85940121
    [[Variables]]
        amplitude:  25.6417887 +/- 0.50569636 (1.97%) (init = 10)
        centerx:    2.00326033 +/- 0.01163397 (0.58%) (init = 2.051478)
        centery:   -0.49692376 +/- 0.01680902 (3.38%) (init = -0.478231)
        sigmax:     0.51074241 +/- 0.01027264 (2.01%) (init = 1)
        sigmay:     1.11741198 +/- 0.02190772 (1.96%) (init = 2)
        rotation:   0.48130689 +/- 0.01542118 (3.20%) (init = 0.1)
    [[Correlations]] (unreported correlations are < 0.100)
        C(centerx, centery)  = +0.5767
        C(amplitude, sigmax) = +0.3342
        C(amplitude, sigmay) = +0.3046




.. GENERATED FROM PYTHON SOURCE LINES 139-140

The process of making plots to check it worked is the same as before.

.. GENERATED FROM PYTHON SOURCE LINES 140-166

.. code-block:: Python

    fig, axs = plt.subplots(2, 2, figsize=(10, 10))

    vmax = np.nanpercentile(Z, 99.9)

    ax = axs[0, 0]
    art = ax.pcolor(X, Y, Z, vmin=0, vmax=vmax, shading='auto')
    plt.colorbar(art, ax=ax, label='z')
    ax.set_title('Data')

    ax = axs[0, 1]
    fit = model.func(X, Y, **result.best_values)
    art = ax.pcolor(X, Y, fit, vmin=0, vmax=vmax, shading='auto')
    plt.colorbar(art, ax=ax, label='z')
    ax.set_title('Fit')

    ax = axs[1, 0]
    fit = model.func(X, Y, **result.best_values)
    art = ax.pcolor(X, Y, Z-fit, vmin=0, vmax=10, shading='auto')
    plt.colorbar(art, ax=ax, label='z')
    ax.set_title('Data - Fit')

    for ax in axs.ravel():
        ax.set_xlabel('x')
        ax.set_ylabel('y')
    axs[1, 1].remove()
    plt.show()



.. image-sg:: /examples/images/sphx_glr_example_two_dimensional_peak_004.png
   :alt: Data, Fit, Data - Fit
   :srcset: /examples/images/sphx_glr_example_two_dimensional_peak_004.png, /examples/images/sphx_glr_example_two_dimensional_peak_004_3_00x.png 3.00x
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 3.144 seconds)


.. _sphx_glr_download_examples_example_two_dimensional_peak.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: example_two_dimensional_peak.ipynb <example_two_dimensional_peak.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: example_two_dimensional_peak.py <example_two_dimensional_peak.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: example_two_dimensional_peak.zip <example_two_dimensional_peak.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
