[gdal-dev] Layer intersection: two polygons with a common point

[Even Rouault]

Le jeudi 28 avril 2016 08:50:46, Ari Jolma a écrit :
> Jukka pointed out this
> 
> http://gis.stackexchange.com/questions/191336/gdal-ogr-ogr-layer-intersecti
> on-producing-points-from-linestrings
> 
> to me.
> 
> It is true that intersecting two layers, where one contains
> 
> 'POLYGON (( 140 360, 140 480, 220 480, 220 360, 140 360 ))'
> 
> and the other
> 
> 'POLYGON (( 220 260, 220 360, 300 360, 300 260, 220 260 ))'
> 
> produces an empty layer. The reason is a check
> 
> (z_geom->IsEmpty() ||
>    (x_geom->getDimension() == 2 &&
>     y_geom->getDimension() == 2 &&
>     z_geom->getDimension() < 2))
> 
> (x and y are the polygons and z is the intersection, a point in this
> case), which leads to discarding the point.
> 
> I believe this check was introduced as a result of ticket 4772
> 
> https://trac.osgeo.org/gdal/ticket/4772
> 
> in changeset 25427, where it is stated that "when intersecting 2
> polygons, do not consider a point or linear intersection as a valid one".

Oh I see we are both involved in this ;-) The advantage of the current 
behaviour is that it is compatible with most popular output formats that don't 
like mixing geometry types in the same layer (think of shapefiles, and 
generally postgis, spatialite when layer geometry type is not unknown, ...). 
Perhaps one possible approach would be to check the output layer geometry type 
: if it is wkbUnknown, keep all non empty geometries, otherwise only keep 
those compatible of the geometry type. And perhaps combine that default 
behaviour with an option to the method "KEEP_LOWER_DIMENSION_GEOMETRIES" 
(would be honoured only if output layer geometry type == wkbUnknown) ? Since I 
feel that in most cases and in real-life, intersection of polygons resulting 
in points or lines should not happen.

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com