[gdal-dev] Question: Layer intersection in Python

[Michaelis, Luis]

Hey Even,

thank you for the exceptionally quick and helpful reply! Switching to the GPKG/MEMORY drivers did indeed fix the problem, and since I don’t rely on the output being a shapefile this works perfectly fine for me. I never thought about switching to a different driver.

Thanks again for the great work!

- Luis

From: Even Rouault <even.rouault at spatialys.com>
Sent: 20 June 2024 19:43
To: Michaelis, Luis <luis.michaelis at iee.fraunhofer.de>; gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] Question: Layer intersection in Python


Luis,

the issue is that shapefile layers do not really support a genuine unknown layer geometry type. They accept it on layer creation, but as soon as you write the first feature into it, the feature layer geometry type is used as the layer geometry type, and shapefile don't accept a mix of points and multipoints. Which brings to another issue is that PROMOTE_TO_MULTI until now only affected lines and polygons, not points. I've fixed this with https://github.com/OSGeo/gdal/pull/10258 .

One workaround for you is to write your result in a Memory layer or a GeoPackage layer for example, and then use gdal.VectorTranslate() to convert it to a multipoint shapefile

Even
Le 20/06/2024 à 14:11, Michaelis, Luis via gdal-dev a écrit :
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




_______________________________________________

gdal-dev mailing list

gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/gdal-dev

--

http://www.spatialys.com

My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20240621/915d5ed7/attachment-0001.htm>