2

I have a shapefile containing multiple polygons of protected zones, I try to compute a raster with, for each pixel, the distance to the nearest protected zone.

I've been able to obtain it with QGIS using the function Distance to nearest hub, with the result looking like (white=protected, red = farthest from protected zone): Proximity to protected areas

I'm now looking for a pythonic way to do this, I didn't find the equivalent function in the Python QGIS API. I did find multiple functions for finding the nearest neighbor for one point, which I could loop over each point, but I beleive there could be better ways to achieve this? What's the best option to scale it to each point?

6

1 Answer 1

1

Thanks for the many interesting methods shared in comments. I decided to call QGIS functions from python directly as I needed other QGIS functions as well.

I had some troubles to have the correct GQIS installation without UI and calling the functions from a script, so I share here how I built my docker image :

from python:3.9.7
RUN apt update
RUN apt-get install -y wget && rm -rf /var/lib/apt/lists/*
RUN wget \
    https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh \
    && mkdir /root/.conda \
    && bash Miniconda3-latest-Linux-x86_64.sh -b \
    && rm -f Miniconda3-latest-Linux-x86_64.sh

ENV PATH="/root/miniconda3/bin:${PATH}"
ARG PATH="/root/miniconda3/bin:${PATH}"
RUN conda update --all -y
RUN conda install -c conda-forge gdal
RUN conda install -c conda-forge qgis==3.16.3
RUN apt-get update -y && apt-get install libgl1 -y
ENV PYTHONPATH="/root/miniconda3/pkgs/qgis-3.16.3-py39h31e92eb_1/share/qgis/python:${PYTHONPATH}"
ENV PYTHONPATH="/root/miniconda3/pkgs/qgis-3.16.3-py39h31e92eb_1/share/qgis/python/plugins:${PYTHONPATH}"

And calling the function like this:

from qgis.core import *
QgsApplication.setPrefixPath('/usr', True)

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

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

grid = QgsVectorLayer('grid.geojson', 'Grid')
hubs = QgsVectorLayer('hubs.geojson', 'Hubs')

parameters = {
    "INPUT": grid,
    "HUBS": joined,
    "FIELD": "WDPA_PID",
    "UNIT": 0,
    "OUTPUT": "/mnt/dev/workspace/protected-areas.geojson"
}
out = Processing.runAlgorithm("qgis:distancetonearesthubpoints", parameters) 

This way I could obtain the same transformations as I did throught QGIS UI!

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.