[gdal-dev] Question: Layer intersection in Python

[Michaelis, Luis]

Hello GDAL devs!
I've got a question about layer intersection detection.

I'd like to use the GDAL Python interface to do batched intersection checks between two vector layers, one containing only two-point line strings and one containing terrain contour lines. To compute these intersections, I would like to use the `ogr.Layer.Intersection` function, since it is supposed do batched intersection detection in native code.

Since all of the intersections will either be points or multi points, I was expecting to use the `KEEP_LOWER_DIMENSION_GEOMETRIES` and `PROMOTE_TO_MULTI` options, to record the intersections into a new layer with a multi-point geometry. This, however, does not work and I simply don't get any intersections at all. Setting the geometry to just points does also not work. In neither of these cases are any errors or warnings issued.

When I use a destination layer with a `wkbUnknown` geometry however, I do get warnings about it trying to add muti-point geometries to a point layer (which I would sort of expect). Setting `SKIP_FAILURES` and saving the resulting shape, I can actually see the intersection points but as expected, I can't see the multipoint intersections. They're simply missing.

So, the question is: am I missing something? Do need to add some metadata information to the destination layer before computing the intersections or is what I'm trying to do impossible?

Here's the code I've written: https://gist.github.com/lmichaelis-fhg/3ae5ba63b92c1a6ad609baa3ccbfbf80#file-dem-generate-rix-gdal-py-L109-L121

- Luis

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20240620/d50d7b69/attachment.htm>