In Python, what is the best way to replicate the QGIS Heatmap plugin?

Question

I would like to create a spatial weighted heatmap in Python where I have control over the boundingbox, grid size and bandwidth. In example if I want to create a population heatmap on a grid of 200*200 meters with a bandwidth of 500 meters:

Using the QGIS Heatmap plugin:

kde(locations = xy, weight = population, boundingbox, gridsize = 200, bandwidth = 500, kernel ="gaussian")

I have not came across a package which are able to do just this.

Following the example: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html From what I understand is that one should first create a meshgrid and the reshape the kde onto the grid. This is what I tried so far, but there are 3 issues:

This part takes very long: zz2 = kernel((np.ravel(x_mesh), np.ravel(y_mesh))) - 470.44 seconds on 19150 points where in Qgis it is mere seconds.
The y-axis is still inverted.
Not sure how to set the bandwidth to a constant 500

Code:

from scipy import stats
from shapely.geometry import Point
import geopandas as gpd
import numpy as np
from osgeo import gdal
from osgeo import osr
import time
import pandas as pd
#Input variables
grid_size=200
h=500
#Data
x =  np.array([-285815.24600105, -285905.88928823, -285596.62853068,
       -285376.49911475, -284530.02007635, -285976.25971212,
       -285079.67702268, -286188.5497945 , -284810.5502149 ,
       -285707.6207015 , -285072.46928953, -284872.60260027,
       -285567.26057971, -284593.23417313, -285318.32010344,
       -285767.26258091, -284600.84807157, -285185.11331713,
       -284727.6299865 , -284982.16195329, -284983.76372273,
       -284753.27862336, -284688.7406417 , -284963.14336973,
       -285102.43887492, -284610.34171822, -284710.3065015 ,
       -284501.4211114 , -286246.95919243, -284923.81296141,
       -285880.14147568, -285099.55526278, -284611.19426662,
       -286108.759291  , -285358.72069313, -284906.19046438,
       -286177.47753297, -284571.34168874, -285519.67954529,
       -285162.43056364, -285915.21656255, -285599.37350284,
       -284494.30220736, -284577.61017269, -284793.26653895,
       -285115.45608425, -285915.80558585])
y = np.array([2906143.2521925 , 2906369.43984717, 2906356.32381486,
       2906525.61255684, 2906540.60531809, 2906586.42258352,
       2906450.18112564, 2906707.0534267 , 2906492.11113259,
       2906725.89380165, 2906763.87804504, 2906779.45155159,
       2906947.06861677, 2906866.70425729, 2906864.30448599,
       2907483.92004085, 2907011.01133657, 2907183.5411114 ,
       2907125.55907197, 2907137.32092455, 2907403.91453819,
       2907417.71324586, 2907309.79221579, 2907636.60663656,
       2907754.1172582 , 2907559.26299843, 2907747.3226264 ,
       2907750.19855555, 2907966.22491989, 2907996.07814694,
       2908213.91807075, 2908003.55703708, 2908039.44317742,
       2908125.88796091, 2908214.69867858, 2908353.7416716 ,
       2908438.61892689, 2908267.34531307, 2908373.17285713,
       2908369.37610769, 2908494.01196971, 2908602.17039364,
       2908507.03090379, 2908737.87072884, 2908685.12160762,
       2908645.49069608, 2908723.21635992])
weight_value = np.array([7985585., 7985084., 7985237., 7984908., 7985446., 7985504.,
       7984242., 7984369., 7984735., 7985019., 7984076., 7984041.,
       7983581., 7984401., 7985564., 7983173., 7984675., 7984697.,
       7984507., 7984368., 7984972., 7984348., 7985082., 7983377.,
       7984336., 7984319., 7984419., 7984460., 7984684., 7984942.,
       7984028., 7985162., 7984346., 7983969., 7984232., 7985258.,
       7984913., 7985284., 7984889., 7984567., 7984341., 7984907.,
       7984793., 7982612., 7983755., 7984752., 7983938.])
#create geodataframe
df_geometry = [Point(xy) for xy in zip(x, y)]
gdf_centroid = gpd.GeoDataFrame(df_geometry, geometry=df_geometry)
gdf_centroid['weigth'] = weight_value
#Create GRID
gdf_centroidg_bb = gdf_centroid.total_bounds
xmin = gdf_centroidg_bb[0]
xmax = gdf_centroidg_bb[2]
ymin = gdf_centroidg_bb[1]
ymax = gdf_centroidg_bb[3]
x_grid = np.arange(xmin-h, xmax+h, grid_size)
y_grid = np.arange(ymin-h, ymax+h, grid_size)
x_mesh, y_mesh = np.meshgrid(x_grid, y_grid)
#Create Kernel Density Estimation
positions = np.vstack([x_mesh.ravel(), y_mesh.ravel()])
values = np.vstack([x, y])
kernel = stats.gaussian_kde(values, weights = weight_value)
kernel.set_bandwidth(bw_method=kernel.factor / 3.)
#This takes too long. (470.44 seconds on 19150 points)
start = time.time()
zz2 = kernel((np.ravel(x_mesh), np.ravel(y_mesh)))
end = time.time()
print(end - start)
#Reshape the kde
zz2 = np.reshape(zz2.T, x_mesh.shape)
#Setup the raster metadata
nrows,ncols = np.shape(y_mesh)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)
#Export kernel density to geotiff
output_raster = gdal.GetDriverByName('GTiff').Create('population_heatmap2.tif',ncols, nrows, 1 ,gdal.GDT_Float32)
output_raster.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(2051)
output_raster.SetProjection( srs.ExportToWkt() )
output_raster.GetRasterBand(1).WriteArray(zz2)
output_raster.FlushCache()

Is there a better why of doing this in a Python script without using QGIS?

score 2 · Answer 1 · answered Sep 20 '20 at 10:17

2

You can call QGIS modules/plugins etc, from outside of QGIS. https://docs.qgis.org/3.4/fi/docs/pyqgis_developer_cookbook/intro.html#using-pyqgis-in-standalone-scripts Also this link talks about using processing algorithms from the console, which may also be of assistance. https://docs.qgis.org/3.10/en/docs/user_manual/processing/console.html#processing-console

Maybe you could try and work out how to call the Heatmap plugin via pyQGIS, essentially replicating the GUI tool outside of QGIS. Check out this https://docs.qgis.org/3.10/en/docs/user_manual/processing_algs/qgis/interpolation.html#python-code

Looks like you can call the algorithm from python, so maybe have a read of the above and see how you go!

answered Sep 20 '20 at 10:17

nr_aus

3,535
6
26

1

Just run the algorithm and in the processing menu there is an option for history. Copy the run command and paste it in the python console. You can modify the run algorithm to read attributes from the table. – kartoza-geek Sep 20 '20 at 20:11
Thanks for the reply. I will explore this option some more, but at the moment I have a stand alone anaconda/python installation from where I run scripts. I don't know if I miss something, but as I understand in order to run qgis.core I need to run it from python that is bundled with qgis. Also, I don't see where I can specify the extent of the heatmap. I want to be able to built different heatmaps where the rasters can be stacked on top of each other. – user19349 Sep 21 '20 at 07:57
There is a facility to run QGIS Python libraries standalone, so no, you dont have to run python from within QGIS.. The link above goes into detail about how to do this (just import the PyQGIS modules as per other python module imports). So it should be possible to integrate this into your existing scripts. – nr_aus Sep 21 '20 at 08:30
Thank you again, I can see how this is supposed to help but I've bee struggling to load qgis.core. I have followed many examples off adding environment variables etc, but keep on getting errors of either qgis module not found or missing dll's. – user19349 Sep 22 '20 at 10:47
1

There could be a few different problems - There are a few articles about this on stack exchange... Try here to begin with https://gis.stackexchange.com/questions/141990/running-python-outside-of-qgis-importerror-no-module-named-qgis-core One of the keys is to ensure consistency with your python version and QGIS libraries ie: 2.7 won't work with QGIS 3.x libraries. – nr_aus Sep 23 '20 at 02:05
1

Thanks nr_aus. I went through many examples trying to resolve the issues stated. The link you provided were the most helpful. Not sure yet how the .pth files should look like exactly, but for now I copied those links to my PYTHONPATH and created a new environment which resolved the qgis.core issue. I will see when I have time to update my question with the process I followed. – user19349 Sep 23 '20 at 11:42

score 1 · Answer 2 · answered Sep 25 '20 at 14:13

I finally managed to run a heatmap using the Qgis modules. For a non-technical user it was a bit of a struggle but definitely worth the effort. Ran a heatmap on points scattered over the whole of Africa on a 200 meter grid in just over 3 minutes.

Software:

Anacondo/Python 3.7.4
Pycharm
Qgis 3.14

Process followed:

To import qgis.core I followed the advice here: enter link description here Basically saying to first add the following to your PATH environment variables:

C:\Program Files\QGIS 3.14\bin;C:\Program Files\QGIS 3.14\apps\qgis\bin

And secondly to create a .pth file in your virtual environment under Lib\site-packages, which contains the lines:

C:\Program Files\QGIS 3.14\apps\qgis\python
C:\Program Files\QGIS 3.14\apps\Python37\
C:\Program Files\QGIS 3.14\apps\Python37\lib\
C:\Program Files\QGIS 3.14\apps\Python37\lib\site-packages
C:\Program Files\QGIS 3.14\bin
C:\Program Files\QGIS 3.14\include
C:\Program Files\QGIS 3.14\apps\qgis\bin

My screenshot of the path:

Next issue was to set the QgsApplication which presented an error: “This application failed to start because it could not find or load the Qt platform plugin” This was resolved by adding the following environment variable to the script:

os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = 'C:\Program Files\QGIS 3.14\apps\Qt5\plugins'

see:

QGIS 3.14 vs "no Qt platform plugin could be initialized".

Now the QgsApplication can be initialised:

from qgis.core import (
     QgsApplication,
     QgsProcessingFeedback,
     QgsVectorLayer
)
QgsApplication.setPrefixPath('C:\Program Files\QGIS 3.14\apps\qgis\', True)
from PyQt5 import QtGui, QtCore
qgs = QgsApplication([], False)
qgs.initQgis()

Next issue was that my processing was missing all native algorithms. To solve this I had to add the following code:

import sys sys.path.append('C:\Program Files\QGIS 3.14\apps\qgis\python\plugins')

import processing from processing.core.Processing import Processing Processing.initialize()

from qgis.analysis import QgsNativeAlgorithms QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

see: > Using QGIS3 Processing algorithms from standalone PyQGIS scripts (outside of GUI)

Note: You may receive an error telling that it cant find the proj.db file. Even with this error the process finished as expected but this may also be resolved by setting the PROJ_LIB environment variable to point to the PROJ.4 data directory (where proj.db lives). see: > Ogr2ogr: ERROR 1: PROJ: pj_obj_create: Cannot find proj.db

This however did not work for me. I added the PROJ_DEBUG = 3 to my environment variables and saw that it searches for the proj.db under C:/Users/user_name/AppData/Roaming/python\profiles\default/proj\proj.db

So I just copied the proj.db to that location.

Finally I could run the heatmap algorithm:

params = {'INPUT':'path to shapefile', 'RADIUS':500, 'RADIUS_FIELD':'', 'PIXEL_SIZE':200, 'WEIGHT_FIELD':'SAL_ID', 'KERNEL':0, 'DECAY':0, 'OUTPUT_VALUE':0, 'OUTPUT':'output\test6.tif'}

processing.run("qgis:heatmapkerneldensityestimation", params)

Full Script:

import os
os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = 'C:\\Program Files\\QGIS 3.14\\apps\\Qt5\\plugins'
#Tried to import the proj.db but did not work
os.environ['GDAL_DATA'] = '/home/server/anaconda3/share/gdal'
os.environ['PROJ_LIB'] = '/home/server/anaconda3/share/proj'
from qgis.core import (
     QgsApplication,
     QgsProcessingFeedback,
     QgsVectorLayer
)
QgsApplication.setPrefixPath('C:\Program Files\QGIS 3.14\apps\qgis\', True)
from PyQt5 import QtGui, QtCore
qgs = QgsApplication([], False)
qgs.initQgis()
import sys
sys.path.append('C:\Program Files\QGIS 3.14\apps\qgis\python\plugins')
import processing
from processing.core.Processing import Processing
Processing.initialize()
from qgis.analysis import QgsNativeAlgorithms
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
params = {'INPUT':'Path to Shape',
          'RADIUS':500,
          'RADIUS_FIELD':'',
          'PIXEL_SIZE':200,
          'WEIGHT_FIELD':'SAL_ID',
          'KERNEL':0,
          'DECAY':0,
          'OUTPUT_VALUE':0,
          'OUTPUT':'output\test6.tif'}
processing.run("qgis:heatmapkerneldensityestimation", params)

In Python, what is the best way to replicate the QGIS Heatmap plugin?

2 Answers2

os.environ['GDAL_DATA'] = '/home/server/anaconda3/share/gdal'

os.environ['PROJ_LIB'] = '/home/server/anaconda3/share/proj'