Thursday, February 15, 2018

Web Mapping Tutorial in QGIS

Web Mapping might be a difficult thing to do for most GIS user, because it needs web programming skills such as Javascript, HTML, PHP, etc. On the other hand web mapping is a nice tool to visualize GIS information and can be shared and accessed easily by many people through a web page. This tutorial will discuss how to make a simple web mapping without any single code in QGIS. For this tutorial I am using QGIS 2.18.

Installing Web Mapping Plugin

The magic behind this web mapping tutorial is QGIS2Web Plugin which will export map layers in QGIS canvas into a web page. Therefore, firstly we have to install this plugin. To install QGIS2Web plugin can be done as follow.
1. In the toolbar menu. Open Plugin > Manage and Install Plugins. As in figure 1.
  
Manage and Install Plugins QGIS
Figure 1. Manage and Install Plugins
2. In the search toolbar, type qgis2web. The plugin in will appear and then check it. Clik Install Plugin at the right bottom of Plugins window.

Install qgis2web plugin
Figure 2. Install qgis2web plugin

Add Data to QGIS Map Canvas

In this tutorial, I used earthquake data that can be downloaded from USGS Earthquake Catalog. Download the data in csv format or geojson. Add the data into QGIS map canvas as in figure 3.

Earthquake point dataset in QGIS
Figure 3. Earthquake point dataset
Open layer properties and do some styling styles. For example I classified the magnitude of earthquake into five classes. So it looks like in figure 5. Furthermore you can make a heatmap or add another data to the map.

QGIS Layer Properties Window
Figure 4. Layer Properties Window
Earthquake point classification based on magnitude
Figure 5. Earthquake point classification based on magnitude

Create A Web Mapping Application

If the map is ready, then it's time to build a web mapping application using qgis2web plugin. The plugin can be found in the web menu, as in figure 6. Soon you will see the plugin window like figure 7.
qgis2web plugin
Figure 6. qgis2web plugin






Export to web map window QGIS
Figure 7. Export to web map window
There are some main sections in this plugin, in particular in the export tab.
  1. Layers and Groups section. It control how a layer will be displayed in the web map, including the pop up fields option. In this part we can also cluster a set points. 
  2. Data Export. This part gives some options for exporting the map into a web map, such as: Location of exported web map files, mapping library location, minify GeoJSON file, etc.
  3. Scale/Zoom.  This section control the map extent including zoom setting.
  4. Appearance. In this section we can add some tools to the web map such as: address searching, layer list, measure tool, pop up, user geolocation. Moreover we can set the extent of the web map in this part.
Export to web map window is our playing field to set some parameters for web mapping application. Must be remember, every time you change or set a parameter, you must click Update Preview button to see the result. Now let's play!

1. Firstly we have to choose web mapping library OpenLayers or Leaflet. For this case I used Leaflet because I think it is lighter and faster.
2. Our web mapping does not have any basemap. Let's pick one by choosing a basemap. For example I chose Open Street Map (OSM) in the list at right bottom. Then clicked Update preview. The web map will show like figure 8. There are some basemap available, just explore them and see the result which is fitted with your purpose.

Web Mapping with OSM basemap QGIS
Figure 8. Web Mapping with OSM basemap
3. In layers and group section. I want to limit some fields to show in the pop up window. Unfortunately we can't do it directly in qgis2web plugin. So go back to QGIS main window and hide the field. You can hide fields in the layer properties. In this case I just want to show time, earthquake magnitude and location/place. Therefore I hide all other fields, except those three one. (To go back to QGIS main window, you don't have to close qgis2web window, just switch it)

Hide a fields in layer properties
Figure 9. Hide a fields in layer properties
4. After done in QGIS main window for hiding some fields. Switch back to qgis2web plugin window. Set the visible fields from no label to inline label, to give each field's name in the pop window. To see the result clik Update preview. The result as in figure 10.

QGIS Web map with pop up window
Figure 10. Web map with pop up window
Other parameters could be set as below.
Data export
Exporter : Export to Folder
Mapping library location: Local
Minify GeoJSON files: Checked (True)
Precision: Maintain

Scale/Zoom
Extent: Fit to layers extent
Max zoom level: 28
Min zoom level: 1
Restric to extent: Checked (True)

Appearance
Add address search: Checked(True)
Add layer list: Expanded
Layer search: place (Field name)
Measure tool: Metric
Template: Fullscreen

In figure 11 can be seen those settings in the plugin window. You can play with them and don't forget to click Update preview to see the result.

Final qgis2web plugin setting for web mapping
Figure 11. Final qgis2web plugin setting for web mapping
If you satisfy with the setting and how the map looks like. Last step click Export button. The exported web map will be opened in a browser as in figure 12 below.

Figure 12. Web Mapping is ready
That's all the tutorial how to create a web mapping application in QGIS. It's quite easy and no coding at all. The web map is working smoothly. See the action in the video below. If you want to explore more about qgis2web plugin, visit the plugin github wiki page here.



Tuesday, January 23, 2018

Creating Heatmap From Scratch in Python

Heatmap is frequently used to visualize event occurrence or density. There are some Python libraries or GIS software/tool that can be used to create a heatmap like QGIS, ArcGIS, Google Table Fusion, etc. Unfortunately, this post won't discussed how to create a heatmap using those software/tool, but more than that, we will write our own code to create a heatmap in Python 3 from scratch using Python common library.
The algorithm which will be used to create a heatmap in Python is Kernel Density Estimation (KDE). Please refer to this post (QGIS Heatmap Using KDE Explained) to get more explanation about KDE and another post (Heatmap Calculation Tutorial) which give an example how to calculate intensity for a point from a reference point using KDE.

Importing Library

Actually, there are some libraries in Python that can be used to create heatmap like Scikit-learn or Seaborn. But we will use just some libraries such as matplotlib, numpy and math. So we are starting with importing those three libraries.

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

Heatmap Dataset

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]

Heatmap Point Dataset in Python
Figure 1. Point Dataset

Grid Size and Radius

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

Getting X,Y Min/Max to Construct Grid

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)

Calculate Grid Center Point

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)

Kernel  Density Estimation Function

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

Compute Density Value for Each Grid

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)

Visualize The Result

The last part we visualize the result using matplotlib color mesh. We also add a color bar to see the intensity value. The heatmap result can be seen in figure 2.

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

Heatmap Output in Python
Figure 2. Heatmap Output
The complete code snippet can be found below.

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()

That's all how to create heatmap in Python from scratch using KDE. There are other kernel shape available like Gaussian, Triweight, Epanechnikov, Triangular, etc. Please free to add those kernel shape and modify the code. Try to experiment with changing some parameters like radius and grid size and explore the result.

Sunday, January 21, 2018

Heatmap Calculation Tutorial Using Kernel Density Estimation (KDE) Algorithm

Kernel Density Estimation (KDE) Overview

The previous post had discussed about Kernel Density Estimation (KDE) in creating a heatmap in QGIS. It explained about background and conceptual approach how KDE is applied for a heatmap production. This post will give a tutorial and example how to calculate a density value estimation around a point dataset. We will go step  by step, so hopefully this tutorial could give you a clear understanding how to do a density calculation using KDE.

As explained before in the previous post, there are some KDE shapes which are using to construct a Probability Density Function (PDF) in order to estimate a density at a location to a reference point. In this tutorial, we will use quartic kernel shape. The quartic kernel shape has a function as in equation 1, with the shape of the kernel can be seen in figure 1. If a density value also considered a weight (W), a constant (K) and Intensity (I), then the function become as in equation 2. In this case we assume W=K=I=1. Therefore we will use equation 1 through out this example.
Quartic Kernel Shape Function
Equation 1. Quartic KDE function

Equation 2. Quartic KDE function with KWI
Quartic Kernel Shape
Figure 1. Quartic kernel shape

Kernel Density Estimation (KDE) Basic Calculation Example

Using the kernel, then we will calculate an estimation density value at a location from a reference point. Figure 2 shows more detail about the quartic kernel shape and some properties such as bandwidth (h), reference point (O), estimation point (z) and the distance (d) from reference point to estimation point. Moreover figure 3 shows how it looks like in a space.

Quadratic Kernel Shape and related properties
Figure 2. Quadratic Kernel Shape and related properties

portrait of KDE properties is space
Figure 3. KDE properties is space
Our goal is to find the density value at z. Knowing the distance (in this example 0.5), then we can calculate the density value using equation 1 as follow.

It is quite simple, isn't it? So we can calculate a density at any point within the bandwidth radius. When a point reaches the bandwidth radius (equal to 1) the density value will be 0, and we do not consider any point outside the bandwidth radius.

Kernel Density Estimation (KDE): More Complex Calculation Example

The example above shows the basic concept in estimating a density at a location from a reference point. Now let's see a more complex example. Figure 4, consist of 3 points data with coordinates O1(6,6)m, O2(10,10)m, O3(5,11)m. Next we want to compute the density estimation value at point z1(7,5)m, z2(7,11)m and z3(7,9)m. The kernel radius is 4 m.
KDE Three points dataset and three estimation points
Figure 4. Three points dataset and three estimation points
Before we do calculation, observe that z2 is in intersection area between O2 and O3 bandwidth radius and z3 is in intersection of all reference points radius. When a point lies in an intersection radius we have to sum up all density value from each reference point, then the density function becomes as in equation 3.
Summation of Quartic KDE function
Equation 3. Summation of Quartic KDE function
 Now, let's do the calculation. Firstly we compute the density at z1, therefore we have to compute the distance O1-z1, with equation 3.
Equation 3. Formula to calculate distance between two points
 Then the distance O1-z1:
 
Great, we just need to use the computed distance, and done! But wait, it can't be like that. The distance can't be used right away. Why? The quartic kernel density function that we're using is in standardize form. It means the bandwidth radius has a fix number as 1. So we have to divide the computed distance with the actual kernel bandwidth (in this case is 4 m). 1.41/4=0.35. Thus the density at z1:



Next we compute the z2. Because z2 is in intersection of O2 and O3, first we calculate the distance of z2 to O2 and O3, then compute the density of z2 from each reference point and sum up the both result. The calculation can be seen as follow:






Lastly for z3:








In this tutorial we just calculated density for some points around the dataset. In a real case, the algorithm will calculate the density for a whole dataset area. For this we have to specify the size of output pixel. The calculation then will take place on each pixel. As the pixel size gets smaller, will produce a smoother result. But on the other hand, will require more resources for computation. Hopefully you could get a robust understanding how a heatmap is created using KDE. You can create a heatmap in QGIS using QGIS heatmap plugin. If you are interested to see how to apply this algorithm in Python see this post (How to Create Heatmap in Python from Scratch)    

Monday, November 20, 2017

QGIS Heatmap Using Kernel Density Estimation Explained

Heatmap is a nice visualization method to display event density or occurrence. Heatmap is also used in clustering points where more points in an area will have higher value compare to less point in the same area. Therefore, with a heatmap we can see a concentration of event's occurrence. For example criminal events in a city can be visualized in a heatmap, so one knows how the crime events spread around the city. Such map is commonly called a crime map.

To create a heatmap is quite easy and straight forward. There are many software and tools that can be  used to produce a heatmap like QGIS, ArcGIS, Crimestat, Google table fusion, etc. We just need to upload point dataset, setting some parameters and the result will come up. But the purpose of this post is not just to produce a heatmap with some clicks. But more than that to explore and trying to get more understanding what is going on behind the screen for data and parameters that have been specified. In another term, how the algorithm  is working in producing a heatmap.

QGIS is an open source GIS software that can be used to produce a heatmap from a set of data point with Heatmap Plugin. The plugin is using Kernel Density Estimation algorithm for creating a heatmap. Because of that I will discuss how this algorithm(Kernel Density Estimation) is applied to process an input point dataset into a heatmap.

Creating A Heatmap in QGIS

To get an insight for what we will discuss, firstly let's create a heatmap of theft from vehicle incident across Vancouver, Canada.

1. Download crime data of Vancouver city 2016 from Vancouver city data catalogue. For convenient I provided the data here. (Actually the data consist of several crime events such as mischief, theft from vehicle, break and enter residential and other theft). I just selected theft from vehicle crime.
2.  Add the data into QGIS map window as in figure 1.

Theft from vehicle accident in Vancouver city 2016
Figure 1. Theft from vehicle accident in Vancouver city 2016
3. Select Heatmap plugin. It can be found in menu Raster >> Heatmap >> Heatmap...(See figure 2)

QGIS Heatmap plugin
Figure 2. Heatmap plugin
4. The heatmap plugin window will appear as in figure 3. Make sure the Input point layer is the theft_fr_vehicle. Give a name for Output raster. I named it theft_fr_vehicle_heatmap. Then set Radius 500. Cliked Advanced option, set Cell size X 20 and Cell size Y 20. Leave the Kernel shape as Quartic(biweight) and the other parameter as default.

QGIS Heatmap Plugin window
Figure 3. Heatmap Plugin window
5. After everything is done. Push OK button and the heatmap will appear as in figure 4.

Theft from vehicle hetmap
Figure 5. Theft from vehicle heatmap
6. Hmm...it looks ugly, doesn't it? So let's make it more colorful. Right click the heatmap layer. In Layer Properties window as in figure 6, choose Style. Change Render type to Singleband pseudo color. Then choose your favorite color in the Color option.

Change heatmap color style
Figure 6. Change heatmap color style
 8. Next goes to figure 7. Still in Layer Properties, here we change the blending mode. Choose Darken (You can experiment with other blending option). Click apply to see the result immediately in the thumbnail preview. Finally push OK button when done. The result of heatmap with the new style can be seen in the figure 8. 

Change Heatmap color blend
Figure 7. Change Heatmap color blend
Theft from vehicle heatmap with singleband pseudo color
Figure 8. Theft from vehicle heatmap with singleband pseudocolor
It's done. The theft from vehicle heatmap of Vancouver city shows us how the incidents are spread in the city. In the north of the city can be seen the color is darker which indicates high density value. So take care when you park a car there :).

Kernel Density Estimation Algorithm

As I mentioned earlier. The heatmap was created with Kernel Density Estimation algorithm. Now let's explore how this algorithm is working, so we can tune related parameters to get a more meaningful heatmap cause we understand how the result comes up.

Kernel Shape

Estimation is predicting an unknown value at a location from reference points. Actually in predicting the unknown value we interpolate the value from known points value. Therefore we also called KDE as Kernel Density Interpolation. In estimation a point value, KDE uses a Probability Density Function (PDF). What is PDF and how could it be used to estimate an unknown value? To answer this question, I will give an illustration. Assume there is a road segment and I want to know how long it is. So I measure it with a tape, first measurement (forward)  I got 100.251 m.Then backward I got 100.310 m? Hmm...why they are different. I am curious and take the third measurement, got the figure 100.279 m. How could this happen? The answer is Error. In every measurement there is an error. Nobody knows a real measurement value except God. But what could I do to get a value closest to the real value? Take more measurement!!! More measurement, more number of  values that close the real value, thus it has more probability to get the real value which is my expected value. But I don't want to waste my time measuring the road 1000 times.......

The same problem already happened hundred years ago. In the 17th century, Galileo made an astronomical observation and noted that the errors from the observation were symmetric and small errors occurred more frequently than large errors. This led to several hypothesized of distribution errors. Later in 1809 Gauss developed the distribution normal function and showed that errors were fitted well by this distribution. (Read the history here). The function of normal distribution, which is also called Gauss distribution can be seen in equation 1. If this function were plotted with some defined parameters, we can see how it looks like in figure 9.

Distribution Normal Equation
Equation 1. Distribution Normal function
Figure 9. Normal distribution curve
The normal distribution curve looks like a bell, and simply in KDE we call this as Kernel Shape. If we look at heatmap plugin, there are some Kernel shapes available, there are: Quartic, Triangular, Uniform and Epanechnikov.  But there are more kernel shapes available like Cosine, Gaussian, Tricube, etc. Some of these shapes can be seen in figure 10.
Kernel Shapes
Figure 10. Kernel Shapes (Source: https://en.wikipedia.org/)
Next question. How those shapes are used in estimation of unknown point? As we already knew, the density value is defined by a Kernel shape with a function . To apply it for a heatmap, the function must be transformed into a "geospatial sense" form, because it will be used to predict an unknown value at a point from a point with known value. Here the distance between known to unknown point is a parameter which must be involved in the function. But where we will put it in the function? For example, in the normal distribution function in equation 1, we can see there is (𝚡-𝜇) term. This term is the distance between a value (x) with the expected value (𝜇). This term can be translated as the distance between known to unknown point (d). Then 𝜎 is standard deviation. Standard deviation in statistics measures variation of dataset values. Actually it is the average distance of x to 𝜇. So this term in KDE is translated into bandwidth (h). For Gaussian, h remains as standard deviation but for other kernels h is radius. Hence the equation of KDE with Gaussian Kernel shape has the form as in equation 2, with the visual illustration can be seen in figure 11. If we standardize the h parameter into 1 (100%), the equation 2 can be simplified into equation 3.
Gaussian KDE function
Equation 2. Gaussian KDE function
Equation 3. Standardize Gaussian KDE Function
Gaussian Kernel Shape
Figure 11. Gaussian Kernel Shape
From figure 11 we can observe that the density value of unknown point will decrease smoothly following the Gaussian Probability Density Function (PDF). It means the unknown value will deviate largely as it moves away from a known point. Other Kernel Shape functions are summarized in the following table.

KernelFunction Shape
Gaussian


Quartic

Epanechnikov

Triangular

Triweight

Uniform

 

How to Choose A "Right" Kernel Shape

This is the hardest part in this post. How to choose a right kernel shape? Honestly I don't know the exact answer. But after some research I found a hint in determining a kernel shape. Based on Crimestat III workbook, there are two questions to answer in deciding a kernel shape. First. How large an area will have the effect from an incidents. Second. How much of the effect will remain at original location and how much will be dispersed throughout the bandwidth interval. For example if you think the effect of accident will decrease in linear pattern from an original location, then Triangular Kernel would be a good option. Otherwise if the accident will have the same effect throughout a specified radius, the Uniform kernel would be a good option. So understanding the two questions above for a case will give a better decision in choosing a kernel shape and related parameters.

The following table gives kernel shape option for some types of crime. Although it is for crime mapping, hopefully it could give the idea in defining a kernel shape for other cases. The table was taken from Crimestat III workbook.

How to choose a kernel shape for estimation

I hope this post could give you an understanding what is a Kernel Density Estimation and its application in creating a heatmap. In the next post, I give a numerical example how this algorithm is used to estimate an unknown value from a reference point. So you will get a complete understanding from theory and its implementation.