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

Thu Apr 28 04:17:52 PDT 2016

Le jeudi 28 avril 2016 12:03:55, Jukka Rahkonen a écrit :
> 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-inters
> > > ecti 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,

OGR layer algebra methods can receive specific options (so strictly speaking 
not general configuration options), and that's already implemented.

e.g for Intersection():
 * The recognized list of options is:
 * <ul>
 * <li>SKIP_FAILURES=YES/NO. Set to YES to go on, even when a
 *     feature could not be inserted or a GEOS call failed.
 * <li>PROMOTE_TO_MULTI=YES/NO. Set to YES to convert Polygons
 *     into MultiPolygons, or LineStrings to MultiLineStrings.
 * <li>INPUT_PREFIX=string. Set a prefix for the field names that
 *     will be created from the fields of the input layer.
 * <li>METHOD_PREFIX=string. Set a prefix for the field names that
 *     will be created from the fields of the method layer.
 * <li>USE_PREPARED_GEOMETRIES=YES/NO. Set to NO to not use prepared
 *     geometries to pretest intersection of features of method layer
 *     with features of this layer.
 * <li>PRETEST_CONTAINMENT=YES/NO. Set to YES to pretest the
 *     containment of features of method layer within the features of
 *     this layer. This will speed up the method significantly in some
 *     cases. Requires that the prepared geometries are in effect.
 * </ul>

> 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

That would be ProperIntersection() then. Why not (though internally we'd want 
to factor code between Intersection and ProperIntersection), but...

GEOS has a GEOSPreparedContainsProperly() to do that, but the concept of 
ProperIntersection() is yet to be defined and not so obvious, especially for 
lines.

Let's say ?:

ProperIntersection(polyA, polyB) := Intersection(polyA, polyB) if it is a 
polygon, or NULL otherwise
ProperIntersection(poly, line) :=
 * Intersection(poly, line) if it is a line (perhaps with the additional 
constraint that the interior of the line must not be contained by the boundary 
of poly, so as to reject lines that are partly or fully in the boundary of 
polygons ? But I'm thinking there are cases where the line would be for 
example made of a segment that is part of the boundary of the poly and another 
segement inside the poly, hence the intersection with the poly would be the 
line itself, but we would want to keep only the part that is in the interior 
of the poly. And I can't think of a way to do that easily with GEOS.) 
  * NULL otherwise ??? (that might be a point)
ProperIntersection(poly, point) := Intersection(poly, point) (perhaps with the 
additional constraint that the point must not be contained by the boundary of 
poly ?) 

For consistency with the above, we should normally have :
ProperIntersection(lineA, lineB) := Intersection(lineA, lineB) if it is a 
line, or NULL otherwise
But use cases for this are not so obvious ( lines that "overlap" in part )

Or perhaps:
ProperIntersection(lineA, lineB) :=
   * Intersection(lineA, lineB) if it is a line,
   * Intersection(lineA, lineB) if it is a point that is contained in the 
interior of both lineA and lineB (ie a point that is not an end of one of  the 
geom)
   *or NULL otherwise

ProperIntersection(line, point) := Intersection(line, point)

So as to generalize the rule would be :

ProperIntersection(geomA, geomB) :=
 * Intersection(geomA, geomB) if dimension of this intersection = min( 
dimension(geomA), dimension(geomB) ), 
 * NULL otherwise

or perhaps:

ProperIntersection(geomA, geomB) :=
  * Intersection(geomA, geomB) if dimension of this intersection = min( 
dimension(geomA), dimension(geomB) ), 
 * Intersection(geomA, geomB) if dimension of this intersection < min( 
dimension(geomA), dimension(geomB) ) but its interior is contained in the 
interior of both geomA and geomB,  ( I'm pretty sure that the only possibility 
with geometries on the 2D plan is when geomA and geomB are lines )
 * NULL otherwise

Hum the more I'm thinking about the rules, the less obvious they are.

I'm still thinking that some option to ensure we don't try to feed the output 
layer with an incompatible geometry would be useful.

> 
> -Jukka Rahkonen-
> 
> 
> 
> 
> 
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev

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