I have a dataset with many overlapping polygons. I need to extract height values from the underlying DEM for each [area] in each polygon, and compare those to properties of the corresponding polygon.

Thus the way forward is to first generate points (pixel centroids) inside polygons, then drape and extract z values, which I use for further comparison.

However, when I generate those pixel centroids, they don't keep any properties from the corresponding polygon. I cannot do a spatial join either, because as stated before, the polygons overlap. Random points don't do the trick either, since their positioning doesn't correspond to the underlying DEM.

Is there some way to join generated points with their corresponding polygon?

Picture for easier understanding (the squares corresponding to background DEM, o = below threshold set by each polygon, x above and to be cut):


  • Do they share the same CRS? Can you share your data with us?
    – Taras
    Commented Jun 20, 2022 at 8:22
  • Cannot really share the data, since the datasets are quite large (>10gb), all data share the same crs Commented Jun 20, 2022 at 8:23
  • Could those points be generated randomly? Beside the "Drape (set Z value from raster)", there is also v.drape
    – Taras
    Commented Jun 20, 2022 at 8:25
  • I wouldn't know how, they need to be regular (as in: one point per pixel of the DEM, for each polygon) Basically each polygon has a threshold value field, and I want to cut out each point with height above the threshold to produce new polygons later on Commented Jun 20, 2022 at 8:30
  • I think there are some mixed concepts here, and it is not clear what you need to achieve. Forget about points and centroids, you need to convert your polygons to raster using the desired threshold value (previously added as a field) making sure to use the same resolution and have the new rasters snapped to your DEM. From here you want to perform a conditional operation in these new rasters and keep the pixels that meet your criteria to convert them back to polygons again
    – Albert
    Commented Jun 20, 2022 at 9:16

1 Answer 1


Found a (dirty) solution without having to resort to iterating over thousands of layers:

  1. Generate points (pixel centroids) inside the polygons
  2. Buffer (10m)
  3. Dissolve by poly_id (generate points creates an own ID list based on the input polygons)
  4. Join attributes by location - layers in [original polygons] which are within [dissolved]
  5. Join by field value, input 1: centroids, poly_id, input 2: spatially joined, poly_id

yielding the point layer with attached attribute table from the original polygons

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.