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@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
