title: Smoothing my heatmap in Python
body: <p>I am developing a script in order to make heatmap from a sky survey with python and the libraries <code>numpy, astropy</code>.</p>
  
  <p>I created a stars distribution map and now I'm trying to make a heatmap.
  My heatmap is done and works well, but my next step is to smooth it with a Gaussian. That's to say, convolute the data by a Gaussian with dispersion = 2 sigma.</p>
  
  <p>The problem is I don't get a good smoothing heatmap. As you could see after, my plots are not good with the convolution function by <code>scipy</code> and/or <code>astropy</code> (I scripted both methods).</p>
  
  <p>This is my code:</p>
  
  <pre><code># -*- coding: utf-8 -*-
  #!/usr/bin/env python
  
  from astropy.io import fits
  from astropy.table import Table
  from astropy.table import Column
  from astropy.convolution import convolve, Gaussian1DKernel
  import numpy as np
  import scipy.ndimage as sp
  import matplotlib.pyplot as plt
  
  
  
          ###################################
          # Importation du fichier de champ #
          ###################################
  
  filename = '/home/valentin/Desktop/Field169_combined_final_roughcal.fits_traite_traiteXY_traiteXY_final'
  
  print 'Fichier en cours de traitement' + str(filename) + '\n'
  
  # Ouverture du fichier à l'aide d'astropy
  field = fits.open(filename)        
  
  # Lecture des données fits 
  tbdata = field[1].data            
  
  
          #######################################
          # Parametres pour la carte de densité #
          #######################################
  
  # Boite des étoiles bleues :
  condition_1 = np.bitwise_and(tbdata['g0-r0'] &gt; -0.5, tbdata['g0-r0'] &lt; 0.8 )  # Ne garder que les -0.4 &lt; (g-r)0 &lt; 0.8
  condition_final = np.bitwise_and(tbdata['g0'] &lt; 23.5, condition_1)       # Récupere les valeurs de 'g0' &lt; 23.5 dans les valeurs de blue_stars_X
  
  Blue_stars = tbdata[condition_final]
  
  RA_Blue_stars = Blue_stars['RA']                        # Récupere les valeurs de 'RA' associées aux étoiles bleues
  DEC_Blue_stars = Blue_stars['DEC']                      # Récupere les valeurs de 'DEC' associées aux étoiles bleues
  
  
  # Boite des étoiles tres bleues :
  condition_2 = np.bitwise_and(tbdata['g0-r0'] &gt; -0.5, tbdata['g0-r0'] &lt; 0.2 )
  condition_final2 = np.bitwise_and(tbdata['g0'] &lt; 23.5, condition_2)
  
  Very_Blue_stars = tbdata[condition_final2]
  
  RA_Very_Blue_stars = Very_Blue_stars['RA']                      # Récupere les valeurs de 'RA' associées aux étoiles bleues
  DEC_Very_Blue_stars = Very_Blue_stars['DEC']
  
  # ==&gt; La table finale avec le masque s'appelle Blue_stars &amp; Very_Blue_stars
  
          ##################################################################
          # Traçage des différents graphiques de la distribution d'étoiles #
          ##################################################################
  
  
  fig1 = plt.subplot(2,2,1)
  plt.plot(tbdata['g0-r0'], tbdata['g0'], 'r.', label=u'Etoiles du champ')
  plt.plot(Blue_stars['g0-r0'], Blue_stars['g0'], 'b.', label =u'Etoiles bleues')
  plt.plot(Very_Blue_stars['g0-r0'], Very_Blue_stars['g0'], 'k.', label =u'Etoiles tres bleues')
  plt.title('Diagramme Couleur-Magnitude')
  plt.xlabel('(g0-r0)')
  plt.ylabel('g0')
  plt.xlim(-1.5,2.5)
  plt.ylim(14,28)
  plt.legend(loc='upper left')
  plt.gca().invert_yaxis()
  
  fig1 = plt.subplot(2,2,2)
  plt.plot(RA_Blue_stars, DEC_Blue_stars, 'b.', label =u'Etoiles bleues', alpha=0.15)
  plt.title('Carte de distribution des etoiles bleues')
  plt.xlabel('RA')
  plt.ylabel('DEC')
  plt.legend(loc='upper left')
  
  fig1 = plt.subplot(2,2,3)
  plt.plot(RA_Very_Blue_stars, DEC_Very_Blue_stars, 'r.', label =u'Etoiles tres bleues',alpha=0.4)
  plt.title('Carte de distribution des etoiles tres bleues')
  plt.xlabel('RA')
  plt.ylabel('DEC')
  plt.legend(loc='upper left')
  
  fig1 = plt.subplot(2,2,4)
  plt.plot(RA_Blue_stars, DEC_Blue_stars, 'b.', label =u'Etoiles bleues', alpha=0.15)
  plt.plot(RA_Very_Blue_stars, DEC_Very_Blue_stars, 'r.', label =u'Etoiles tres bleues',alpha=0.4)
  plt.title('Carte de distribution des etoiles bleues et tres bleues')
  plt.xlabel('RA')
  plt.ylabel('DEC')
  plt.legend(loc='upper left')
  
          ######################################################################
          # Traçage des différents graphiques de la carte de densité d'étoiles #
          ######################################################################
  
  ###############################################################################
  # Carte de densité des étoiles bleues pour 1 pixel de 1 arcmin^2 (bins = 180) #
  ###############################################################################
  
  X_Blue_stars = Blue_stars['X']
  Y_Blue_stars = Blue_stars['Y']
  
  heatmap, xedges, yedges = np.histogram2d(X_Blue_stars, Y_Blue_stars, bins=180) # bins de 180 car 3° de champ en RA = 180 arcmin de champ en RA
  extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
  
  plt.clf()
  plt.subplot(2,2,1)
  plt.imshow(heatmap, extent=extent)
  plt.colorbar()
  plt.title('Carte de densite des etoiles bleues (non lisse)')
  plt.xlabel("X")
  plt.ylabel("Y")
  plt.gca().invert_xaxis()
  
  
  ####################################################################################################################################
  # Carte de densité lissée (par convolution avec une gaussienne 2 sigma) des étoiles bleues pour 1 pixel de 1 arcmin^2 (bins = 180) #
  # ==&gt; Avec Scipy                                                        #
  ####################################################################################################################################
  
  lissage_X_scipy = sp.filters.gaussian_filter(X_Blue_stars, sigma = 2, order = 0)
  lissage_Y_scipy = sp.filters.gaussian_filter(Y_Blue_stars, sigma = 2, order = 0)
  
  
  heatmap_lisse_scipy, xedges_lisse_scipy, yedges_lisse_scipy = np.histogram2d(lissage_X_scipy, lissage_Y_scipy, bins=180)
  extent_lisse_scipy = [xedges_lisse_scipy[0], xedges_lisse_scipy[-1], yedges_lisse_scipy[0], yedges_lisse_scipy[-1]]
  
  plt.subplot(2,2,2)
  plt.imshow(heatmap_lisse_scipy, extent=extent_lisse_scipy)
  plt.colorbar()
  plt.title('Carte de densite des etoiles bleues lisse a 2 sigma (scipy)')
  plt.xlabel("X")
  plt.ylabel("Y")
  plt.gca().invert_xaxis()
  
  ####################################################################################################################################
  # Carte de densité lissée (par convolution avec une gaussienne 2 sigma) des étoiles bleues pour 1 pixel de 1 arcmin^2 (bins = 180) #
  # ==&gt; Avec Astropy                                                          #
  ####################################################################################################################################
  
  # Creation du kernel :
  
  K = Gaussian1DKernel(stddev=2) # Détermination de la déviation standard (sigma)
  
  lissage_X_astropy = convolve(X_Blue_stars, kernel=K, boundary='fill')
  lissage_Y_astropy = convolve(Y_Blue_stars, kernel=K, boundary='fill')
  
  heatmap_lisse_astropy, xedges_lisse_astropy, yedges_lisse_astropy = np.histogram2d(lissage_X_astropy, lissage_Y_astropy, bins=180)
  extent_lisse_astropy = [xedges_lisse_astropy[0], xedges_lisse_astropy[-1], yedges_lisse_astropy[0], yedges_lisse_astropy[-1]]
  
  plt.subplot(2,2,3)
  plt.imshow(heatmap_lisse_astropy, extent=extent_lisse_astropy)
  plt.colorbar()
  plt.title('Carte de densite des etoiles bleues lisse (astropy)')
  plt.xlabel("X")
  plt.ylabel("Y")
  plt.gca().invert_xaxis()
  plt.show()
  
  print "Création du Diagramme"
  </code></pre>
  
  <p>I get this :</p>
  
  <p><img src="https://i.stack.imgur.com/nC2Db.jpg" alt=""></p>
  
  <ol>
  <li>Top left: heatmap without smoothing</li>
  <li>Top right: heatmap with scipy smoothing</li>
  <li>Bottom: heatmap with astropy smoothing</li>
  </ol>
  
  <p>I don't know why, there are lots of holes, lacks ... with the smoothing :/</p>
  
  <p><strong>UPDATE:</strong> </p>
  
  <p>After the answer from Framester, I wrote an easier script which contains the "same thing" that my problem. I applied the same method (by scipy for example) and I get a smoothing heatmap :)</p>
  
  <pre><code>import matplotlib.pyplot as plt
  import numpy as np
  import scipy.ndimage as sp
  
  
  x = np.random.randn(100000)
  y = np.random.randn(100000) + 5
  
  # normal distribution center at x=0 and y=5
  fig1 = plt.subplot(2,2,1)
  plt.hist2d(x, y, bins=40)
  plt.colorbar()
  plt.title('Heatmap without smoothing')
  plt.xlabel("X")
  plt.ylabel("Y")
  
  
  # smoothing
  
  X = sp.filters.gaussian_filter(x, sigma = 2, order = 0)
  Y = sp.filters.gaussian_filter(y, sigma = 2, order = 0)
  
  
  heatmap, xedges, yedges = np.histogram2d(X, Y, bins=40)
  extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
  
  fig1 = plt.subplot(2,2,2)
  plt.imshow(heatmap, extent=extent)
  plt.colorbar()
  plt.title('Heatmap with smoothing')
  plt.xlabel("X")
  plt.ylabel("Y")
  plt.show()
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/w9W8p.png" alt="Not smooth and smooth"></p>
  
  <p>So, my question is : why my script doesn't work ? :/</p>
  
  <p><strong>SOLUTION BY MSEIFERT :</strong></p>
  
  <pre><code>plt.clf()
  plt.subplot(2,2,1)
  plt.imshow(heatmap, extent=extent, interpolation='none')
  plt.colorbar()
  plt.title('Carte de densite des etoiles bleues (non lisse)')
  plt.xlabel("X")
  plt.ylabel("Y")
  plt.gca().invert_xaxis()
  
  plt.subplot(2,2,2)
  plt.imshow(convolve(heatmap, Gaussian2DKernel(stddev=2)), interpolation='none')
  plt.colorbar()
  plt.title('Carte de densite des etoiles bleues lisse (astropy)')
  plt.xlabel("X")
  plt.ylabel("Y")
  plt.gca().invert_xaxis()
  
  plt.show()
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/4OsGH.png" alt="Final result"></p>