2

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:

enter image description here

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:

  1. 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.
  2. The y-axis is still inverted.
  3. 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?

0
1

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!

6
  • 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. 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 7: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 8:30
  • 1
    There could be a few different problems - There are a few articles about this on stack exchange... Try here to begin with gis.stackexchange.com/questions/141990/… 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 2: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
0

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:

  1. 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: enter image description here


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

  1. 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.


  1. 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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.