#### Importing Library

import matplotlib.pyplot as plt import numpy as np import math

To create a heatmap, we need a point dataset that consist of x,y coordinates. Here we create two list for x and y. The plot of dataset can be seen in figure 1.

#POINT DATASET x=[20,28,15,20,18,25,15,18,18,20,25,30,25,22,30,22,38,40,38,30,22,20,35,33] y=[20,14,15,20,15,20,32,33,45,50,20,20,20,25,30,38,20,28,33,50,48,40,30,35]

Figure 1. Point Dataset |

In creating heatmap using KDE we need to specify the bandwidth or radius of the kernel shape and output grid size. For this case, I use radius 10 m and grid size 1 m. Later you can change these parameters to see how they affect the heatmap result.

#DEFINE GRID SIZE AND RADIUS(h) grid_size=1 h=10

To construct grid we use mesh grid. Therefore we need to find x,y minimum and maximum to generate a sequence number of x and y. These sequence numbers then will be used to construct mesh grid. To include all the dataset coverage with a little bit more space, I subtract x,y minimum with radius and add it up for x,y maximum.

#GETTING X,Y MIN AND MAX x_min=min(x) x_max=max(x) y_min=min(y) y_max=max(y) #CONSTRUCT GRID x_grid=np.arange(x_min-h,x_max+h,grid_size) y_grid=np.arange(y_min-h,y_max+h,grid_size) x_mesh,y_mesh=np.meshgrid(x_grid,y_grid)

After constructing mesh grid. Next we calculate the center point for each grid. This can be done with adding x mesh and y mesh coordinate with half of grid size. The center point will be used later to calculate the distance of each grid to dataset points.

#GRID CENTER POINT xc=x_mesh+(grid_size/2) yc=y_mesh+(grid_size/2)

To calculate a point density or intensity we use a function called kde_quartic. We are using Quartic kernel shape, that's why it has "quartic" term in the function name. This function has two arguments: point distance(d) and kernel radius (h).

#FUNCTION TO CALCULATE INTENSITY WITH QUARTIC KERNEL def kde_quartic(d,h): dn=d/h P=(15/16)*(1-dn**2)**2 return P

This is the hardest part of this post. Computing the density value for each grid. We are doing this in three looping. First loop is for mesh data list or grid. Second loop for each center point of those grids and third loop to calculate the distance of the center point to each dataset point. Using the distance, then we compute the density value of each grid with kde_quartic function which already defined before. It will return a density value for each distance to a data point. Here we only consider the point with a distance within the kernel radius. We do not consider the point outside the kernel radius and set the density value to 0. Then we sum up all density value for a grid to get the total density value for the respective grid The total density value then is stored in a list which is called intensity_list.

#PROCESSING intensity_list=[] for j in range(len(xc)): intensity_row=[] for k in range(len(xc[0])): kde_value_list=[] for i in range(len(x)): #CALCULATE DISTANCE d=math.sqrt((xc[j][k]-x[i])**2+(yc[j][k]-y[i])**2) if d<=h: p=kde_quartic(d,h) else: p=0 kde_value_list.append(p) #SUM ALL INTENSITY VALUE p_total=sum(kde_value_list) intensity_row.append(p_total) intensity_list.append(intensity_row)

#HEATMAP OUTPUT intensity=np.array(intensity_list) plt.pcolormesh(x_mesh,y_mesh,intensity) plt.plot(x,y,'ro') plt.colorbar() plt.show()

Figure 2. Heatmap Output |

import matplotlib.pyplot as plt import numpy as np import math #POINT DATASET x=[20,28,15,20,18,25,15,18,18,20,25,30,25,22,30,22,38,40,38,30,22,20,35,33,35] y=[20,14,15,20,15,20,32,33,45,50,20,20,20,25,30,38,20,28,33,50,48,40,30,35,36] #DEFINE GRID SIZE AND RADIUS(h) grid_size=1 h=10 #GETTING X,Y MIN AND MAX x_min=min(x) x_max=max(x) y_min=min(y) y_max=max(y) #CONSTRUCT GRID x_grid=np.arange(x_min-h,x_max+h,grid_size) y_grid=np.arange(y_min-h,y_max+h,grid_size) x_mesh,y_mesh=np.meshgrid(x_grid,y_grid) #GRID CENTER POINT xc=x_mesh+(grid_size/2) yc=y_mesh+(grid_size/2) #FUNCTION TO CALCULATE INTENSITY WITH QUARTIC KERNEL def kde_quartic(d,h): dn=d/h P=(15/16)*(1-dn**2)**2 return P #PROCESSING intensity_list=[] for j in range(len(xc)): intensity_row=[] for k in range(len(xc[0])): kde_value_list=[] for i in range(len(x)): #CALCULATE DISTANCE d=math.sqrt((xc[j][k]-x[i])**2+(yc[j][k]-y[i])**2) if d<=h: p=kde_quartic(d,h) else: p=0 kde_value_list.append(p) #SUM ALL INTENSITY VALUE p_total=sum(kde_value_list) intensity_row.append(p_total) intensity_list.append(intensity_row) #HEATMAP OUTPUT intensity=np.array(intensity_list) plt.pcolormesh(x_mesh,y_mesh,intensity) plt.plot(x,y,'ro') plt.colorbar() plt.show()