
.. 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

        Click :ref:`here <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:: default

    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-27

.. code-block:: default

    npoints = 10000
    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 28-29

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

.. GENERATED FROM PYTHON SOURCE LINES 29-39

.. code-block:: default

    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
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 40-41

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

.. GENERATED FROM PYTHON SOURCE LINES 41-45

.. code-block:: default

    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

 Out:

 .. code-block:: none

    [[Fit Statistics]]
        # fitting method   = leastsq
        # function evals   = 74
        # data points      = 10000
        # variables        = 5
        chi-square         = 23221.7077
        reduced chi-square = 2.32333244
        Akaike info crit   = 8435.02425
        Bayesian info crit = 8471.07595
    [[Variables]]
        amplitude:  27.6616344 +/- 0.68573297 (2.48%) (init = 19.87323)
        centerx:    1.99810718 +/- 0.01471383 (0.74%) (init = 2.123473)
        centery:   -0.49717812 +/- 0.02011906 (4.05%) (init = -0.6082653)
        sigmax:     0.55187429 +/- 0.01262250 (2.29%) (init = 1.66637)
        sigmay:     0.73663336 +/- 0.01696361 (2.30%) (init = 0.8328972)
        fwhmx:      1.29956462 +/- 0.02972371 (2.29%) == '2.3548200*sigmax'
        fwhmy:      1.73463897 +/- 0.03994626 (2.30%) == '2.3548200*sigmay'
        height:     106.882391 +/- 3.63826501 (3.40%) == '1.5707963*amplitude/(max(1e-15, sigmax)*max(1e-15, sigmay))'
    [[Correlations]] (unreported correlations are < 0.100)
        C(amplitude, sigmay) =  0.253
        C(amplitude, sigmax) =  0.220




.. GENERATED FROM PYTHON SOURCE LINES 46-48

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 48-76

.. code-block:: default

    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
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 77-83

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 83-99

.. code-block:: default

    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 100-101

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

.. GENERATED FROM PYTHON SOURCE LINES 101-120

.. code-block:: default

    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
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 121-127

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 127-137

.. code-block:: default

    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

 Out:

 .. code-block:: none

    [[Fit Statistics]]
        # fitting method   = leastsq
        # function evals   = 73
        # data points      = 10000
        # variables        = 6
        chi-square         = 11554.1100
        reduced chi-square = 1.15610466
        Akaike info crit   = 1456.56121
        Bayesian info crit = 1499.82325
    [[Variables]]
        amplitude:  25.4520940 +/- 0.50893773 (2.00%) (init = 10)
        centerx:    1.99532738 +/- 0.01162425 (0.58%) (init = 2.077652)
        centery:   -0.50221822 +/- 0.01657490 (3.30%) (init = -0.5359784)
        sigmax:     0.50576994 +/- 0.01047067 (2.07%) (init = 1)
        sigmay:     1.11968508 +/- 0.02191517 (1.96%) (init = 2)
        rotation:   0.48387508 +/- 0.01534530 (3.17%) (init = 0.1)
    [[Correlations]] (unreported correlations are < 0.100)
        C(centerx, centery)  =  0.537
        C(amplitude, sigmax) =  0.371
        C(amplitude, sigmay) =  0.250




.. GENERATED FROM PYTHON SOURCE LINES 138-139

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

.. GENERATED FROM PYTHON SOURCE LINES 139-165

.. code-block:: default

    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
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_examples_example_two_dimensional_peak.py:


.. only :: html

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



  .. 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-jupyter

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


.. only:: html

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

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