[gdal-dev] Question: Layer intersection in Python

[Even Rouault]

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
> 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/20240620/7bc60a72/attachment.htm>