# Plotting density map of points?

I have some points and I would like to plot the density map of them. I have tried like this:

``````plt.figure()
xedges = np.linspace(540000, 550000, 20)
yedges = np.linspace(5200000, 5500000, 20)
H, xedges, yedges = np.histogram2d(x,y, bins=(xedges, yedges))
print(H)
H = H.T
X, Y = np.meshgrid(xedges, yedges)
plt.pcolormesh(X, Y, H)
plt.colorbar()
`````` ``````from scipy.stats.kde import gaussian_kde
x = np.array(x)
y = np.array(y)
k = gaussian_kde(np.vstack([x, y]))
xi, yi =
np.mgrid[x.min():x.max():x.size**0.5*1j,y.min():y.max():y.size**0.5*1j]
zi = k(np.vstack([xi.flatten(), yi.flatten()]))
fig = plt.figure(figsize=(7,8))
ax1.pcolormesh(xi, yi, zi.reshape(xi.shape))
ax1.set_xlim(x.min(), x.max())
ax1.set_ylim(y.min(), y.max())
ax2.set_xlim(x.min(), x.max())
ax2.set_ylim(y.min(), y.max())
`````` but I don't want my density map would be blocky. I want something smoother like below. If anybody knows what I should do, please let me know. • Both your code examples use coarse and few raster cells, if you want a higher resolution, you need to change that. That's (potentially) all. – bugmenot123 Sep 27 at 14:30

From my understanding, what you looking for is a Kernel Density. Such as KernelDensity in sklearn.neighbors package. You can found the explicit code to create exact example of density map you shown. You can look here : scikit-learn Density Kernel

Code: The following code work perfectly.

``````# Author: Jake Vanderplas <jakevdp@cs.washington.edu>
#

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_species_distributions
from sklearn.datasets.species_distributions import construct_grids
from sklearn.neighbors import KernelDensity

# if basemap is available, we'll use it.
# otherwise, we'll improvise later...
try:
from mpl_toolkits.basemap import Basemap
basemap = True
except ImportError:
basemap = False

# Get matrices/arrays of species IDs and locations
data = fetch_species_distributions()
species_names = ['Bradypus Variegatus', 'Microryzomys Minutus']

Xtrain = np.vstack([data['train']['dd lat'],
data['train']['dd long']]).T
ytrain = np.array([d.decode('ascii').startswith('micro')
for d in data['train']['species']], dtype='int')
Xtrain *= np.pi / 180.  # Convert lat/long to radians

# Set up the data grid for the contour plot
xgrid, ygrid = construct_grids(data)
X, Y = np.meshgrid(xgrid[::5], ygrid[::5][::-1])
land_reference = data.coverages[::5, ::5]

xy = np.vstack([Y.ravel(), X.ravel()]).T
xy *= np.pi / 180.

# Plot map of South America with distributions of each species
fig = plt.figure()

for i in range(2):
plt.subplot(1, 2, i + 1)

# construct a kernel density estimate of the distribution
print(" - computing KDE in spherical coordinates")
kde = KernelDensity(bandwidth=0.04, metric='haversine',
kernel='gaussian', algorithm='ball_tree')
kde.fit(Xtrain[ytrain == i])

# evaluate only on the land: -9999 indicates ocean
Z = Z.reshape(X.shape)

# plot contours of the density
levels = np.linspace(0, Z.max(), 25)
plt.contourf(X, Y, Z, levels=levels, cmap=plt.cm.Reds)

if basemap:
print(" - plot coastlines using basemap")
m = Basemap(projection='cyl', llcrnrlat=Y.min(),
urcrnrlat=Y.max(), llcrnrlon=X.min(),
urcrnrlon=X.max(), resolution='c')
m.drawcoastlines()
m.drawcountries()
else:
print(" - plot coastlines from coverage")
plt.contour(X, Y, land_reference,
levels=[-9998], colors="k",
linestyles="solid")
plt.xticks([])
plt.yticks([])

plt.title(species_names[i])

plt.show()
``````

The result it gives: Also, for an Esri solution, you can follow instruction here ESRI density kernel which I prefer less.