5

Dataset 1: around 80 vectorlayers with polygons. Each polygon has a value for a noise measurement. Some of the layers (up to ten at at time) are intersecting at some locations.

Dataset 2: a vectorlayer with sever hundred points.

The goal: I want to add a column to the point dataset with the maximum value at this place from any of the polygons from the first dataset, as in: "look at the point in question, then look if any of the other layers have a polygon there. Then look at all the values from the noise column of all these polygons and write the highest one into a new column in the attribute set of the point layer".

The problem: If I just merge the noise-layers I later can't use join_attributes_by_location with the point layer as I get a fix geometry error (I guess because the features are overlapping). Do I have to split all the features which overlap with a feature in a different layer by all the features with which they overlap and then merge all layers and then merge all features with the same extent preserving only the maximum noise value for each shape, and then use that layer for join_attributes_by_location? This sounds really complex and also as if it will take forever.

How can I achieve the situation described above (under goal)? Is there another spatial join I can use?

Note 1: I use QGIS 3 and know some Python and very basic R. Note 2: I don't want to find the maximum value in a column in one attribute table, just to clarify.

Note 2: I'm new to SE, this is my first question. I already searched for similar problems but couldn't find anything. Forgive me if I overlooked something or forgot some piece of information.

4

You could use Point Sampling Tool plugin.

Note) If you run Point Sampling Tool plugin on a polygon layer which has overlaps, the tool take the attribute value only from the top feature found at the location.

(1) Merge all polygon layers (which produces a polygon layer with a lot of overlapping polygons).

(2) Run Order by expression tool (Processing Toolbox > Vector general)

enter image description here

Select your noise value field as the expression and tick Sort ascending checkbox, so that the polygon with maximum noise comes on top.

(3) Use Point sampling tool to take the noise value from the newly created Ordered layer.

enter image description here

(4) Open the attribute table of the output from the Point sampling tool.

enter image description here

(... my second sampling point missed all polygons, unfortunately).

  • It didn't quite work, as the geometry was invalid. Also it took forever. But I think your approach is the way to go for smaller datasets than mine. Thank you for the extensive answer! – C4X Nov 28 '19 at 12:50
2

Since you tagged pyqgis, here is a python solution which you can run from the python console. A couple of notes on the script. Since you say have around 80 polygon layers, the script below takes all polygon layers in the current project and stores them in a list, so when you run this script you should make sure your project does not contain any polygon layers which are not part of this 'noise measurement' dataset. Otherwise, if you are handy with python, feel free to devise any other way to access the relevant polygon layers and store them in a list.

You will also need to change the name of the point layer in the code to match yours, as well as the name of the field in your polygon layers which contains the noise measurement data (I hope the field name is the same in every polygon layer, otherwise we're stuffed!). I have commented the code where to change the names.

A new column will be created in the point layer called 'Max_noise' containing the maximum noise reading from all the overlapping polygons. If no polygon geometry intersects the point, the field should be NULL for that feature.

I tested it on a small number of layers with test data and it should do what you want. Hopefully it will not be too slow on your data sets- several hundred points x 80 layers x ?? number of polygons in each layer is quite a few iterations so it might take a while to run.

_project = QgsProject().instance()
pnt_lyr = _project.mapLayersByName('Noise_points')[0] #Change point layer name to match yours
pnt_lyr.dataProvider().addAttributes([QgsField('Max_noise', QVariant.Int)])
pnt_lyr.updateFields()
poly_lyrs = [l for l in QgsProject().instance().mapLayers().values() if isinstance(l, QgsVectorLayer) and l.geometryType() == 2]
polys = []
for l in poly_lyrs:
    polys.append([f for f in l.getFeatures()])
max_fld = pnt_lyr.dataProvider().fields().lookupField('Max_noise')
point_fts = [f for f in pnt_lyr.getFeatures()]
pt_index = QgsSpatialIndex()
for ft in point_fts:
    pt_index.insertFeature(ft)
for f in point_fts:
    id = f.id()
    pt_geom = f.geometry()
    noise_vals = []
    for l in polys:
        for f in l:
            if pt_index.intersects(f.geometry().boundingBox()):
                if f.geometry().intersects(pt_geom):
                    noise_vals.append(f.attribute('Noise_level')) #change field name to match yours
    if noise_vals:
        atts = {max_fld: max(noise_vals)}
        pnt_lyr.dataProvider().changeAttributeValues({id: atts})

Result after running the code. The canvas shows sample points labelled in red with their feature ids, and overlapping polygons labelled with 'Noise_level'. In the point layer's attribute table you can see that the 'Max_noise' field has been created and filled with the maximum value of all overlapping polygons which intersect each point:

enter image description here

1

I would try this approach:

  • Merge all the layers into one (without fusing the polygons)
  • Run the "Union" command to properly split overlapping polygons, while retaining all information.
  • Create a second layer, let's call it "MaxP", with polygons split at each intersection, and only keeping the maximum value found among the original polygons at that position.
  • Reorder the polygons by increasing values (as @Kazuhito suggested, the Point Sample Tool gets the value from the top polygon it encounters, so if you are looking for Max values, use the ascending order).
  • Sample the "MaxP" layer with your points.

Merge the polygon layers into one

Make sure all polygon layers have the same table structure, then go on and merge them into one.

In the Toolbox search "Merge vector layers"

enter image description here

Union

This splits the overlapping polygons into smaller ones. Each overlapping part is duplicated once for each overlapping polygon, so it retains the original information.

enter image description here

Join attributes by location (summarise)

  1. Always in the toolbox, search "Join attributes by location (summarise)" and run it with the following options.
  2. Name of the "Union" layer.
  3. Same: name of the "Union" layer.
  4. Select all geometric predicates of your interest.
  5. Select the name of the field you want to sample (for example "Noise_Level").
  6. Under summaries to calculate, select "max".

enter image description here

Order by expression

This will create the final layer, which we called "MaxP", which you can sample with your points layer. Always in the Processing Toolbox, search "Order by expression" and order the polygons by ascending noise level value.

H/T @Kazuhito

0

So, in the end I solved it using Python, namely Pandas, Geopandas and Shapely. The interesting part is the sorting, duplicated and the spatial join (sjoin) towards the end. The comments in the code should explain most of what is happening. Note: This script creates a table with as many "Max_X"-columns as there are noise-shapefiles. I couldn't figure out how to make a column in the end with a Max value of all the Max_X-values. loc and iloc didn't allow me to slice until the last column like so: loc[3,:]. However, that's probably for a different question.

It's basically doing the same as @Kazuhito's answer suggested, but without splitting the polygons first.

Thank you everyone for the extensive answers! This was my first post on this site after passively reading for a long time and I definitely feel welcome :) Here's the code:

import os

import geopandas as gpd
from geopandas.tools import sjoin
import pandas as pd
from shapely.geometry import Point


# load the subjects measurements as csv
csv = pd.read_csv(
    filepath_or_buffer="/path/to/file.csv",
    usecols=[
        "ID_geo", "A_id", "A_date", "location_accuracy",
        "location_latitude", "location_longitude"
    ],
)

# convert them to a dataframe
df = pd.DataFrame(csv)

# sort by subject iD (A_id), datetime (A_date) and location accuracy (location_accuracy)
df = df.sort_values(["A_id", "A_date", "location_accuracy"])

# create a new column ("duplicated") with boolean values. All rows where the combination
# of A_id and A_date is not unique get "True" unless they are they first occurence of
# this specific tuple. Together with the sorting this can be translated as: "For every
# A_id and A_date, keep only the one value with the lowest location_accuracy (lower is better)"
df["duplicated"] = df.duplicated(subset=["A_id","A_date"], keep="first")

# keep only the ones with False
df = df.query('duplicated == False')

# convert the df to a geodataframe
points = gpd.GeoDataFrame(
    df.drop(['location_latitude', 'location_longitude', 'duplicated'],
    axis=1),
    crs={'init': 'epsg:4326'},
    geometry=[Point(xy) for xy in zip(df.location_longitude, df.location_latitude)]
)

# reproject to EPSG 3035 (the crs of the noise measurement shapefiles)
points = points.to_crs({'init': 'epsg:3035'})

## load and match all shapefiles of noise from Berlin, NRW, Bavaria and state roads
path_to_data = "/path/to/shapefile/folder/"
counter = 0

for filename in os.listdir(path_to_data):
    if filename.endswith(".shp"):
        print(filename)
        # creates a new "DB_Max_X"-str every time the loop runs
        db_max_str = "DB_Max" + str(counter)
        # read the shapefile
        poly = gpd.read_file(os.path.join(path_to_data, filename))
        # set the CRS (set != reproject!)
        poly.crs = {'init': 'epsg:3035'}
        # set the new column to the max value of both columns.
        # in theory, "DB_High" should be higher than "DB_Low"
        # but it isn't reliably
        poly[db_max_str] = poly[["DB_High", "DB_Low"]].max(axis=1)
        # drop the other columns if they exist
        poly_cleaned = poly.drop(
            ['DB_High', 'DB_Low', 'CTRYID', 'Agglomerat', 'UnAggID', 'ICAO'],
            axis=1, 
            errors='ignore'
        )
        # for all the points check if they're within any of the polygons and if so
        # get the values from the polygons
        points = sjoin(points, poly_cleaned, how="left", op="within").drop(
            labels=['index_right', 'index_left'],
            axis=1,
            errors='ignore',
        )
        counter += 1

points.to_csv("/path/to/matched/data.csv")

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.