Trong Python, cách tốt nhất để sao chép plugin QGIS Heatmap là gì?

Question

Tôi muốn tạo một bản đồ nhiệt có trọng số không gian bằng Python trong đó tôi có quyền kiểm soát hộp giới hạn, kích thước lưới và băng thông. Ví dụ, nếu tôi muốn tạo bản đồ nhiệt dân số trên lưới 200 * 200 mét với băng thông 500 mét:

Sử dụng plugin QGIS Heatmap:

kde (vị trí = xy, trọng số = dân số, hộp giới hạn, kích thước lưới = 200, băng thông = 500, nhân = "gaussian")

Tôi chưa tìm thấy một gói nào có thể làm được điều này.

Sau ví dụ: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.htmlTừ những gì tôi hiểu là trước tiên người ta nên tạo một lưới lưới và định hình lại kde trên lưới. Đây là những gì tôi đã thử cho đến nay, nhưng có 3 vấn đề:

1. Phần này mất rất nhiều thời gian: zz2 = kernel ((np.ravel (x_mesh), np.ravel (y_mesh))) - 470,44 giây trên 19150 điểm trong đó Qg là chỉ vài giây.
2. Trục y vẫn ngược.
3. Không chắc chắn cách đặt băng thông thành 500 không đổi

Mã:

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ó lý do nào tốt hơn để thực hiện điều này trong tập lệnh Python mà không sử dụng QGIS không?

Answer 1

Bạn có thể gọi các mô-đun / plugin QGIS, v.v., từ bên ngoài QGIS. https://docs.qgis.org/3.4/fi/docs/pyqgis_developer_cookbook/intro.html#using-pyqgis-in-standalone-scripts Ngoài ra, liên kết này nói về việc sử dụng các thuật toán xử lý từ bảng điều khiển, cũng có thể có ích. https://docs.qgis.org/3.10/en/docs/user_manual/processing/console.html#processing-console

Có lẽ bạn có thể thử và tìm ra cách gọi plugin Heatmap qua pyQGIS, về cơ bản là sao chép công cụ GUI bên ngoài QGIS. Kiểm tra cái nàyhttps://docs.qgis.org/3.10/en/docs/user_manual/processing_algs/qgis/interpolation.html#python-code

Có vẻ như bạn có thể gọi thuật toán từ python, vì vậy có thể đọc phần trên và xem cách bạn thực hiện!

Answer 2

Cuối cùng tôi đã quản lý để chạy một bản đồ nhiệt bằng cách sử dụng các mô-đun Qgis. Đối với một người dùng không chuyên về kỹ thuật, đó là một chút khó khăn nhưng chắc chắn đáng để nỗ lực. Lập bản đồ nhiệt về các điểm nằm rải rác trên toàn châu Phi trên lưới 200 mét chỉ trong hơn 3 phút.

Phần mềm:

• Anacondo / Python 3.7.4
• Pycharm
• Qgis 3,14

Quy trình tiếp theo:

1. Để nhập qgis.core, tôi đã làm theo lời khuyên ở đây: nhập mô tả liên kết ở đây Về cơ bản nói rằng trước tiên hãy thêm phần sau vào các biến môi trường PATH của bạn:

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

Và thứ hai để tạo tệp .pth trong môi trường ảo của bạn trong Lib \ site-package, chứa các dòng:

• 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-package
• C: \ Program Files \ QGIS 3.14 \ bin
• C: \ Program Files \ QGIS 3.14 \ include
• C: \ Program Files \ QGIS 3.14 \ apps \ qgis \ bin

Ảnh chụp màn hình của tôi về đường dẫn:

---

1. Vấn đề tiếp theo là đặt QgsApplication đã xuất hiện lỗi: “Ứng dụng này không thể khởi động vì không thể tìm thấy hoặc tải plugin nền tảng Qt” Điều này đã được giải quyết bằng cách thêm biến môi trường sau vào tập lệnh:

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

xem:

QGIS 3.14 so với "không có plugin nền tảng Qt nào có thể được khởi tạo" .

Bây giờ QgsApplication có thể được khởi tạo:

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. Vấn đề tiếp theo là quá trình xử lý của tôi thiếu tất cả các thuật toán gốc. Để giải quyết vấn đề này, tôi đã phải thêm mã sau:

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

   nhập xử lý từ process.core.Processing import Processing Processing.initialize ()

   từ qgis.analysis nhập QgsNativeAlgorithm QgsApplication.processingRegistry (). addProvider (QgsNativeAlgorithm ())

xem:> Sử dụng thuật toán xử lý QGIS3 từ các tập lệnh PyQGIS độc lập (bên ngoài GUI)

Lưu ý: Bạn có thể nhận được lỗi thông báo rằng nó không thể tìm thấy tệp proj.db. Ngay cả với lỗi này, quá trình đã kết thúc như mong đợi nhưng điều này cũng có thể được giải quyết bằng cách đặt biến môi trường PROJ_LIB trỏ đến thư mục dữ liệu PROJ.4 (nơi proj.db sống). xem:> Ogr2ogr: LỖI 1: PROJ: pj_obj_create: Không thể tìm thấy proj.db

Tuy nhiên, điều này đã không làm việc cho tôi. Tôi đã thêm PROJ_DEBUG = 3 vào các biến môi trường của mình và thấy rằng nó tìm kiếm proj.db trong C: / Users / user_name / AppData / Roaming / python \ profiles \ default / proj \ proj.db

Vì vậy, tôi chỉ sao chép proj.db vào vị trí đó.

---

1. Cuối cùng tôi có thể chạy thuật toán bản đồ nhiệt:

   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'}

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