I'm trying to compare the values of each pixel in relation to each band to create a scatter plot of density. First I extracted the values from the image to form an array, then I created a dataset. Building a pairplot with seaborn I found no problem, now I would like to use sns.Pairgrid in a range from 0 to 255 the result is very poor.This is a test on a clip of the image:

with rasterio.open('C:\data_mosaic.tif', 'r') as ds:
     Arr = ds.read() 
(88, 105, 8)
    b1  b2  b3  b4  b5  b6  b7  b8
0   94.122932   87.930649   57.192429   35.844883   57.971062   65.494003    52.297470  52.553162
1   92.231049   87.693893   53.804562   33.005547   52.124733   56.096642    48.072334  46.176899
2   89.846649   87.448158   49.858879   29.900572   46.716476   44.890785    44.026333  40.420742
3   87.181229   87.291374   46.363262   27.649641   41.478992   36.512981    40.489635  35.537495
4   85.915497   87.230659   43.459812   25.325624   37.368202   30.755083    37.228760  31.470888 
... ... ... ... ... ... ... ... ...

(9240 rows × 8 columns)

g = sns.PairGrid(df)

enter image description here

How do I get to this?

enter image description here

  • Are you trying to do that scatterplot for each subplot in the PairGrid (excluding the diagonal)? – Marcelo Villa-Piñeros Dec 17 '19 at 8:49
  • Yess!!! but I find it very difficult to replicate that diagram for my case...Can you help me? – vins_26 Dec 17 '19 at 19:05

Here is a good answer on how to create a density scatter plot using a Kernel Density Estimate (KDE). Building on top of that answer, and slightly modifying some of the code for readability, here is a function that accomplish the same thing:

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
from matplotlib.colors import Normalize
from scipy.stats import gaussian_kde

def density_scatter_plot(x, y, **kwargs):
    :param x: data positions on the x axis
    :param y: data positions on the y axis
    :return: matplotlib.collections.PathCollection object
    # Kernel Density Estimate (KDE)
    values = np.vstack((x, y))
    kernel = gaussian_kde(values)
    kde = kernel.evaluate(values)

    # create array with colors for each data point
    norm = Normalize(vmin=kde.min(), vmax=kde.max())
    colors = cm.ScalarMappable(norm=norm, cmap='viridis').to_rgba(kde)

    # override original color argument
    kwargs['color'] = colors

    return plt.scatter(x, y, **kwargs)

How should you use this function to create the density scatter plot in each subplot of your PairGrid (excluding the diagonal)? Well, seaborn's PairGrid has a map_offdiag() method which allows you to pass a custom function that accepts both x and y, as well as keyword arguments (**kwargs). This means that it will plot on all the subplots of the grid (again, exlcuding the diagonal) whatever you define in your function. Therefore, your code could look something like this:

g = sns.PairGrid(df)
g = g.map_diag(sns.distplot)
g = g.map_offdiag(density_scatter_plot, s=0.1)  # s controls the size of each point in the scatter plot
g.set(xlim=(0, 255))
g.set(ylim=(0, 255))
g.fig.set_size_inches(10, 10)

Using a random but somewhat similar dataset, I got the following result:

enter image description here

Note that by using **kwargs on the density_scatter_plot() function you can pass whatever keyword argument you want to the plt.scatter() function (e.g. the size s of each point). The map_offdiag() seems to pass a default color argument and that is the reason this argument is overriden inside the density_scatter_plot() function.

  • Thank you for your valuable answer, I would also like to ask you if looking at the dataframe I reported in the question, what is the best solution to define (x - y) considered by your answer, and if there is a way to print only the point cloud rather than the entire background. – vins_26 Dec 18 '19 at 17:02
  • 1
    @vins_26 x and y are passed automatically to the function and correspond to each series in your DataFrame in a given moment (for example b1 and b2 in the second subplot of the first row). The whole point cloud is being plotted (all values from one band vs another band in a given moment). I guess you can define a specific overlap threshold so you just plot the overlapping points but that seems like a different question than the one you asked. – Marcelo Villa-Piñeros Dec 18 '19 at 22:23
  • Thank you so much for your help, your answer helped me!!!! – vins_26 Dec 20 '19 at 19:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.