In Python, qual è il modo migliore per replicare il plugin QGIS Heatmap?

Question

Vorrei creare una mappa termica ponderata spaziale in Python in cui ho il controllo sul riquadro di delimitazione, le dimensioni della griglia e la larghezza di banda. Ad esempio, se voglio creare una mappa di calore della popolazione su una griglia di 200 * 200 metri con una larghezza di banda di 500 metri:

Utilizzando il plugin QGIS Heatmap:

kde (posizioni = xy, peso = popolazione, boundingbox, gridsize = 200, bandwidth = 500, kernel = "gaussian")

Non ho trovato un pacchetto in grado di fare proprio questo.

Seguendo l'esempio: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.htmlDa quello che ho capito è che si dovrebbe prima creare una meshgrid e rimodellare il kde sulla griglia. Questo è quello che ho provato finora, ma ci sono 3 problemi:

1. Questa parte richiede molto tempo: zz2 = kernel ((np.ravel (x_mesh), np.ravel (y_mesh))) - 470,44 secondi su 19150 punti dove in Qgis sono pochi secondi.
2. L'asse y è ancora invertito.
3. Non sei sicuro di come impostare la larghezza di banda su un valore costante di 500

Codice:

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

C'è un motivo migliore per farlo in uno script Python senza usare QGIS?

Answer 1

Puoi chiamare moduli / plugin di QGIS ecc. Dall'esterno di QGIS. https://docs.qgis.org/3.4/fi/docs/pyqgis_developer_cookbook/intro.html#using-pyqgis-in-standalone-scripts Anche questo link parla dell'utilizzo di algoritmi di elaborazione dalla console, che possono anche essere di aiuto. https://docs.qgis.org/3.10/en/docs/user_manual/processing/console.html#processing-console

Forse potresti provare a capire come chiamare il plugin Heatmap tramite pyQGIS, essenzialmente replicando lo strumento GUI al di fuori di QGIS. Controlla questohttps://docs.qgis.org/3.10/en/docs/user_manual/processing_algs/qgis/interpolation.html#python-code

Sembra che tu possa chiamare l'algoritmo da Python, quindi forse leggi quanto sopra e guarda come vai!

Answer 2

Alla fine sono riuscito a eseguire una mappa di calore utilizzando i moduli Qgis. Per un utente non tecnico è stata un po 'difficile, ma ne è valsa la pena. Ho eseguito una mappa di calore su punti sparsi in tutta l'Africa su una griglia di 200 metri in poco più di 3 minuti.

Software:

• Anacondo / Python 3.7.4
• Pycharm
• Qgis 3.14

Processo seguito:

1. Per importare qgis.core ho seguito il consiglio qui: inserisci la descrizione del link qui In pratica dicendo di aggiungere prima quanto segue alle tue variabili d'ambiente PATH:

C: \ Programmi \ QGIS 3.14 \ bin; C: \ Programmi \ QGIS 3.14 \ apps \ qgis \ bin

E in secondo luogo per creare un file .pth nel tuo ambiente virtuale in Lib \ site-packages, che contiene le righe:

• C: \ Programmi \ QGIS 3.14 \ apps \ qgis \ python
• C: \ Programmi \ QGIS 3.14 \ apps \ Python37 \
• C: \ Programmi \ QGIS 3.14 \ apps \ Python37 \ lib \
• C: \ Programmi \ QGIS 3.14 \ apps \ Python37 \ lib \ site-packages
• C: \ Programmi \ QGIS 3.14 \ bin
• C: \ Programmi \ QGIS 3.14 \ include
• C: \ Programmi \ QGIS 3.14 \ apps \ qgis \ bin

Il mio screenshot del percorso:

---

1. Il problema successivo era impostare QgsApplication che presentava un errore: "Questa applicazione non si è avviata perché non è riuscita a trovare o caricare il plugin della piattaforma Qt" Questo è stato risolto aggiungendo la seguente variabile d'ambiente allo script:

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

vedere:

QGIS 3.14 vs "nessun plugin della piattaforma Qt potrebbe essere inizializzato" .

Ora QgsApplication può essere inizializzato:

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. Il problema successivo era che nella mia elaborazione mancavano tutti gli algoritmi nativi. Per risolvere questo problema ho dovuto aggiungere il seguente codice:

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

   importazione elaborazione da processing.core.Processing importazione elaborazione Processing.initialize ()

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

vedere:> Utilizzo degli algoritmi di elaborazione QGIS3 da script PyQGIS standalone (al di fuori della GUI)

Nota: potresti ricevere un messaggio di errore indicante che non è possibile trovare il file proj.db. Anche con questo errore il processo è terminato come previsto, ma questo può anche essere risolto impostando la variabile d'ambiente PROJ_LIB in modo che punti alla directory dei dati PROJ.4 (dove risiede proj.db). vedi:> Ogr2ogr: ERRORE 1: PROJ: pj_obj_create: Impossibile trovare proj.db

Questo però non ha funzionato per me. Ho aggiunto PROJ_DEBUG = 3 alle mie variabili d'ambiente e ho visto che cerca proj.db in C: / Users / nome_utente / AppData / Roaming / python \ profiles \ default / proj \ proj.db

Quindi ho appena copiato il file proj.db in quella posizione.

---

1. Finalmente ho potuto eseguire l'algoritmo della mappa di calore:

   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)

---

Script completo:

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)