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

[Jukka Rahkonen]

Even Rouault <even.rouault <at> spatialys.com> writes:

> 
> 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.
> 


Hi,

I would say that the current behaviour is OK but it should not be called as
"Intersection" because it is so clearly against the OGC Simple Feature
Access - SQL option. There is even a conformity test number 47 on page 63
that takes a linestring and a polygon as inputs and requires a point as an
output. PostGIS, Spatialite and Oracle are also all returning point for this
query.

The test in the standard is:

T47 Intersection(g1 Geometry, g2Geometry) : Geometry
For this test, we will determine the intersection between 
Cam Stream and Blue Lake.
The answer must be 'POINT( 52 18 )'

Instead of adding a new weird configuration options, how about having an OGC
compliant function "Intersects" and our our extension "IntersectsProperly"
like there is ContainsProperly in PosGIS
http://postgis.net/docs/ST_ContainsProperly.html

-Jukka Rahkonen-