# Visualization

This page provides examples of how to visualize membrane curvature via Matplotlib. Two different approaches are suggested.

1. imshow

2. contourf

Warning

Please note that when plotting via `imshow()` or `contourf()`, the orientation of the plot in the final rendering is determined by the `origin` argument. By default, in these Matplotlib functions, the origin is set to the `(left, top)` corner. Additionally, for an array of shape `(M, N)`, the first index `M` runs along the vertical axis, and the second index `N` runs along the horizontal axis.

For these reasons, the returned arrays of surface, mean, and Gaussian curvature should be transposed, and `origin=lower` should be passed to `imshow()` or `contourf()` to generate the correct plots. This is a critical step when visualizing results obtained from MembraneCurvature analyses to avoid generating plots of membrane surface and curvature with the wrong orientation.

## 1. imshow

A simple plot using `imshow()` can be obtained like so:

```import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.imshow(avg_mean_curvature.T, cmap='bwr', interpolation='gaussian', origin='lower')
ax.set_title('Mean Curvature')
plt.show()
```

With `imshow()`, each element of the array is plotted as a square in a matrix of M x N elements. Since we are interested in plotting the array wwhere the first index runs along the horizontal axis, while the second index runs along the vertical, we should transpose the array obtained from the MembraneCurvature analysis. In the generated plot, the color of each square is determined by the value of the corresponding array element and the colormap of preference.

For example, to visualize the results obtained in 1. Membrane-only systems, we can run:

```In [1]: import MDAnalysis as mda
...: from membrane_curvature.base import MembraneCurvature
...: from MDAnalysis.tests.datafiles import Martini_membrane_gro
...: import matplotlib.pyplot as plt
...:

In [2]: universe = mda.Universe(Martini_membrane_gro)

In [3]: curvature_upper_leaflet = MembraneCurvature(universe,
...:                                             select='resid 1-225 and name PO4',
...:                                             n_x_bins=8,
...:                                             n_y_bins=8,
...:                                             wrap=True).run()
...:

In [4]: mean_upper_leaflet = curvature_upper_leaflet.results.average_mean

In [5]: curvature_lower_leaflet = MembraneCurvature(universe,
...:                                             select='resid 226-450 and name PO4',
...:                                             n_x_bins=8,
...:                                             n_y_bins=8,
...:                                             wrap=True).run()
...:

In [6]: mean_lower_leaflet = curvature_lower_leaflet.results.average_mean

In [7]: leaflets = ['Lower', 'Upper']

In [8]: curvatures = [mean_lower_leaflet, mean_upper_leaflet]

In [9]: fig, [ax1, ax2] = plt.subplots(ncols=2, figsize=(6,3), dpi=200)
...: for ax, mc, lf in zip((ax1, ax2), curvatures, leaflets):
...:     ax.imshow(mc.T, origin='lower', interpolation='gaussian', cmap='seismic')
...:     ax.set_aspect('equal')
...:     ax.set_title('{} Leaflet'.format(lf))
...:     ax.axis('off')
...:
```

## 2. contourf

You can use contour plots using `contourf()`. With this approach, contour lines and filled contours of the obtained two-dimensional data are plotted. A contour line connects points with the same curvature values.

When plotting using `contourf()`, an extra step is required to perform an interpolation. We suggest using `scipy.ndimage.gaussian_filter()` as in:

```In [10]: from scipy import ndimage

In [11]: leaflets = ['Lower', 'Upper']

In [12]: fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(5,3))
....: for ax, mc, lf in zip((ax1, ax2), curvatures, leaflets):
....:     arr_ = ndimage.gaussian_filter(mc, sigma=1, order=0, mode='reflect')
....:     ax.contourf(arr_.T,
....:                 cmap='bwr',
....:                 levels=30)
....:     ax.set_aspect('equal')
....:     ax.set_title('{} Leaflet'.format(lf))
....:     ax.axis('off')
....:
```