3

I have a set of points (meteorological stations) with 5km buffers around them. Additionally, I have a set of polygons (municipalities) in which those stations are located.

I would like to have data on how much of the buffer is filled by each polygon but I was not able to find a function or a way to do it. So for example for point A I would have that in the buffer of 5km there is 80% of polygon B and 20% of polygon C as the output.

Edit: The anwsers provided do anwser this question in a setting in which I know what polygons and buffers I need to have in the code or when I can do it manualy. Problem is that in my setting with 1 shapefile having 500 polygons and one vector layer having 13000 buffers I am unable to do this manually and I am asking if there is an easier way to do it. I tried implementing the approaches in the questions but I was unable to use them on my example in big scale. Moreover I had problems with coding that I also noted in comments down below.

Main difference to the other similar question provided is that I have 2 vector layers (Buffer with 500 buffers and Municipalities a shape file with 500 different polygons). Previous questions were analyzing this question in situations in which I know precisely the polygons and buffers analyzed, at least that is how I understood them, so after trying to use those approaches on my datasets I was not able to replicate them.

Another way of approaching this problem for me is to use join attributes by location -one-to-many and then getting averages of the duplicate results, only problem is that then I take average of polygons instead of their weights corresponding to % of intersection with the buffer. In my case I do not know what polygons of the shapefile are going to intersect what buffers, or even if the buffer is going to be intersected by 1,2 or more polygons

  • 2
    Have you tried intersection from processing? – Zoltan Apr 13 at 13:53
  • 2
    Did you try the Vector Analysis - Overlay Analysis tool yet? Alternately you could run a Vector overlay, (e.g. Union) , recalculate the area field and then export the table to a spreadsheet to analyze things. – johns Apr 13 at 13:57
3

In this case, you can do it manually

for example, you have this:

enter image description here

You can calculate the area of the buffer and save in the field area

enter image description here

Now, you calculate the intersect between buffer and polygons, and use the next expression to calculate the percentage:

$area/area * 100

Finally, you have the percentage of the buffer is contained in the polygon.

enter image description here

| improve this answer | |
  • Hey, I have been trying to do what you suggested, but it does not work when I have multiple buffers and municipalities interesecting them, and doing it for each one of them does not work, because there is 13k of them. Sometimes the buffer is within the municipality but sometimes it intersects 3 of them. – voorten Apr 14 at 11:46
2

You can achieve this using the DB Manager of QGIS.

You choose Database / Database Manager / Database Manager then Virtual Layers / Qgis Layers you can then try the following query :

select buf.id, mun.name, st_area(st_intersection(buf.geometry, mun.geometry))/st_area(buf.geometry)
from buffer_layer buf, municipalities_layer mun
where st_intersects(buf.geometry, mun.geometry)
order by buf.id

I assume your buffer table has an 'id' attribute and your municipalities layer a 'name' attribute.

You can then load the result in the layer manager.

| improve this answer | |
  • Thanks for help! I tried this but I have never used DB Manager or SQL. So adapting your code I have select SamplingPointLocalId, JPT_KOD_JE, st_area(st_intersection(Buffered.geometry, Fixed geometries powiaty.geometry))/st_area(Buffered.geometry) from Buffered buf, Fixed geometries powiaty mun where st_intersects(buf.geometry, mun.geometry) order by SamplingPointLocalId But it does not work. I am not really sure what should I put instead of "buf.geometry", "mun.geometry" and mun/buf in "from" line. – voorten Apr 14 at 11:16
2

You can use shapely to intersect the geometries and get the area percentage for each, like this:

enter image description here

import json
from pathlib import Path
import matplotlib.pyplot as plt
from shapely.geometry import shape

ROOT_PATH = Path(__file__).parent / ".data"

AREA_1 = json.load(open(ROOT_PATH / "area1.geojson", "r"))
AREA_2 = json.load(open(ROOT_PATH / "area2.geojson", "r"))
INTERSECTED_AREA = json.load(open(ROOT_PATH / "intersected_area.geojson", "r"))

if __name__ == '__main__':
    geom_1 = shape(AREA_1["features"][0]["geometry"])
    geom_2 = shape(AREA_2["features"][0]["geometry"])
    geom_to_intersect = shape(INTERSECTED_AREA["features"][0]["geometry"])

    total_area = geom_to_intersect.area

    intersected_geom_1 = geom_to_intersect.intersection(geom_1)
    intersected_geom_2 = geom_to_intersect.intersection(geom_2)

    total_area_1 = intersected_geom_1.area
    total_area_2 = intersected_geom_2.area

    total_area_1_percentage = (total_area_1 * 100) / total_area
    total_area_2_percentage = (total_area_2 * 100) / total_area

    cm = plt.get_cmap('RdBu')
    fig, axs = plt.subplots()
    axs.set_aspect('equal', 'datalim')

    xs, ys = geom_1.exterior.xy
    axs.fill(xs, ys, alpha=0.5, fc='b', ec='none')

    xs, ys = geom_2.exterior.xy
    axs.fill(xs, ys, alpha=0.5, fc='r', ec='none')

    xs, ys = geom_to_intersect.exterior.xy
    axs.fill(xs, ys, alpha=0.5, fc='g', ec='none')

    plt.show()

    print(total_area)
    print(total_area_1)
    print(total_area_2)
    print(total_area_1_percentage)
    print(total_area_2_percentage)

The matplot code is just for visualization

| improve this answer | |

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.