1.. _rfc-64:
2
3=======================================================================================
4RFC 64: Triangle, Polyhedral surface and TIN
5=======================================================================================
6
7Authors: Avyav Kumar Singh, Even Rouault
8
9Contact: avyavkumar at gmail dot com, even.rouault at spatialys.com
10
11Status: Adopted, implemented
12
13Implementation version: GDAL 2.2
14
15Summary
16-------
17
18As of now, the :cpp:class:`OGRGeometry` class (the base class from which
19all the subtypes are derived) is limited to OGRCompoundCurve,
20OGRCircularString, OGRLinearRing, OGRMultiLineString, OGRMultiPoint,
21OGRMultiPolygon, OGRMultiCurve, OGRSimpleCurve, OGRCurvePolygon and
22OGRPolygon.
23
24This RFC addresses the addition of the following new geometries in
25OGRGeometry:
26
27-  Triangle - A subset of polygons, the fundamental difference is that
28   it is made of 3 nodes only (actually 4, with the last one being the
29   repetition of the first one) and ONLY ONE exterior boundary and NO
30   interior polygons.
31-  PolyhedralSurface - A 3D figure made exclusively of Polygons.
32-  TriangulatedSurface - A subset of PolyhedralSurface; a 3D figure
33   which consists exclusively of Triangles.
34
35Reference documents
36-------------------
37
38-  `OpenGIS Simple Feature Access Part 1 : Common Architecture,v
39   1.2.1 <http://portal.opengeospatial.org/files/?artifact_id=25355>`__,
40   a.k.a. SFA 1.2.1
41-  `BNF of WKT
42   encoding <https://github.com/postgis/postgis/blob/svn-trunk/doc/bnf-wkt.txt>`__:
43   extracted from SQL/MM Part 3
44-  `BNF of WKB
45   encoding <https://github.com/postgis/postgis/blob/svn-trunk/doc/bnf-wkb.txt>`__:
46   extracted from SQL/MM Part 3
47
48Core changes
49------------
50
51The new class hierarchy is the following and is mostly consistent with
52SQL/MM Part 3
53
54.. image:: ../../../images/rfc64/classOGRGeometry_RFC64.png
55
56Some prelimenary work had already been done prior to this proposal, such
57as including the necessary WKB codes in <ogr_core.h>.
58
59Additionally, the `SFCGAL <http://www.sfcgal.org/>`__ library is a new
60optional dependency of GDAL (build support only done for Unix for now).
61The minimum version tested to build is 1.2.2 (as found in Ubuntu 16.04).
62As mentioned in its home page, "SFCGAL is a C++ wrapper library around
63CGAL with the aim of supporting ISO 19107:2013 and OGC Simple Features
64Access 1.2 for 3D operations." It is mostly used as a potential geometry
65backend by PostGIS. It has a C API, that is the one we use.
66
67SFCGAL functions may be used by methods of OGRGeometry (currently
68IsValid(), Distance(), ConvexHull(), Intersection(), Union(),
69Difference(), SymDifference(), Crosses()), as soon as one of the
70geometry operands is a Triangle, PolyhedralSurface or TIN.
71
72Two new OGRGeometry methods are used to convert SFCGAL geometries <->
73OGR geometries.
74
75::
76
77   static sfcgal_geometry_t* OGRexportToSFCGAL(OGRGeometry *poGeom);
78   static OGRGeometry* SFCGALexportToOGR(sfcgal_geometry_t* _geometry);
79
80Besides SFCGAL, GEOS methods are still used in some cases, but with the
81following limitations - a Triangle is converted to a Polygon with one
82exterior ring; Polyhedral Surfaces and Triangulated Surfaces are
83converted to geometry collection of polygons. (each Triangle in a
84Triangulated Surface is converted to a Polygon as described previously)
85
86The API for the new geometries introduced includes -
87
88-  Overwriting existing methods for Polygon in the case of Triangle API.
89   A complete API is provided below -
90
91::
92
93   class CPL_DLL OGRTriangle : public OGRPolygon
94   {
95     private:
96       bool quickValidityCheck() const;
97
98     protected:
99   //! @cond Doxygen_Suppress
100       virtual OGRSurfaceCasterToPolygon   GetCasterToPolygon() const CPL_OVERRIDE;
101       virtual OGRErr importFromWKTListOnly( char ** ppszInput, int bHasZ, int bHasM,
102                                          OGRRawPoint*& paoPoints, int& nMaxPoints,
103                                          double*& padfZ ) CPL_OVERRIDE;
104   //! @endcond
105
106     public:
107       OGRTriangle();
108       OGRTriangle(const OGRPoint &p, const OGRPoint &q, const OGRPoint &r);
109       OGRTriangle(const OGRTriangle &other);
110       OGRTriangle(const OGRPolygon &other, OGRErr &eErr);
111       OGRTriangle& operator=(const OGRTriangle& other);
112       virtual ~OGRTriangle();
113       virtual const char *getGeometryName() const CPL_OVERRIDE;
114       virtual OGRwkbGeometryType getGeometryType() const CPL_OVERRIDE;
115
116       // IWks Interface
117       virtual OGRErr importFromWkb( unsigned char *, int = -1,
118                                     OGRwkbVariant=wkbVariantOldOgc ) CPL_OVERRIDE;
119       virtual OGRErr importFromWkt( char ** ) CPL_OVERRIDE;
120
121       // New methods rewritten from OGRPolygon/OGRCurvePolygon/OGRGeometry
122       virtual OGRErr addRingDirectly( OGRCurve * poNewRing ) CPL_OVERRIDE;
123
124   //! @cond Doxygen_Suppress
125       static OGRGeometry* CastToPolygon(OGRGeometry* poGeom);
126   //! @endcond
127   };
128
129-  The PolyhedralSurface API is derived from OGRSurface. Internally, it
130   uses an OGRMultiPolygon to store all the Polygons comprising the
131   Polyhedral Surface. Most of the implementations of the methods just
132   reference corresponding OGRMultiPolygon methods with checks to ensure
133   that conditions are maintained.
134
135::
136
137   class CPL_DLL OGRPolyhedralSurface : public OGRSurface
138   {
139     protected:
140   //! @cond Doxygen_Suppress
141       friend class OGRTriangulatedSurface;
142       OGRMultiPolygon oMP;
143       virtual OGRSurfaceCasterToPolygon      GetCasterToPolygon() const CPL_OVERRIDE;
144       virtual OGRSurfaceCasterToCurvePolygon GetCasterToCurvePolygon() const CPL_OVERRIDE;
145       virtual OGRBoolean         isCompatibleSubType( OGRwkbGeometryType ) const;
146       virtual const char*        getSubGeometryName() const;
147       virtual OGRwkbGeometryType getSubGeometryType() const;
148       OGRErr exportToWktInternal (char ** ppszDstText, OGRwkbVariant eWkbVariant, const char* pszSkipPrefix ) const;
149
150       virtual OGRPolyhedralSurfaceCastToMultiPolygon GetCasterToMultiPolygon() const;
151       static OGRMultiPolygon* CastToMultiPolygonImpl(OGRPolyhedralSurface* poPS);
152   //! @endcond
153
154     public:
155       OGRPolyhedralSurface();
156       OGRPolyhedralSurface(const OGRPolyhedralSurface &poGeom);
157       virtual ~OGRPolyhedralSurface();
158       OGRPolyhedralSurface& operator=(const OGRPolyhedralSurface& other);
159
160       // IWks Interface
161       virtual int WkbSize() const CPL_OVERRIDE;
162       virtual const char *getGeometryName() const CPL_OVERRIDE;
163       virtual OGRwkbGeometryType getGeometryType() const  CPL_OVERRIDE;
164       virtual OGRErr importFromWkb( unsigned char *, int=-1, OGRwkbVariant=wkbVariantOldOgc ) CPL_OVERRIDE;
165       virtual OGRErr exportToWkb( OGRwkbByteOrder, unsigned char *, OGRwkbVariant=wkbVariantOldOgc ) const CPL_OVERRIDE;
166       virtual OGRErr importFromWkt( char ** )  CPL_OVERRIDE;
167       virtual OGRErr exportToWkt( char ** ppszDstText, OGRwkbVariant=wkbVariantOldOgc ) const  CPL_OVERRIDE;
168
169       // IGeometry methods
170       virtual int getDimension() const  CPL_OVERRIDE;
171
172       virtual void empty()  CPL_OVERRIDE;
173
174       virtual OGRGeometry *clone() const  CPL_OVERRIDE;
175       virtual void getEnvelope(OGREnvelope * psEnvelope) const  CPL_OVERRIDE;
176       virtual void getEnvelope(OGREnvelope3D * psEnvelope) const  CPL_OVERRIDE;
177
178       virtual void flattenTo2D() CPL_OVERRIDE;
179       virtual OGRErr transform(OGRCoordinateTransformation*) CPL_OVERRIDE;
180       virtual OGRBoolean Equals(OGRGeometry*) const CPL_OVERRIDE;
181       virtual double get_Area() const CPL_OVERRIDE;
182       virtual OGRErr PointOnSurface(OGRPoint*) const CPL_OVERRIDE;
183
184       static OGRMultiPolygon* CastToMultiPolygon(OGRPolyhedralSurface* poPS);
185       virtual OGRBoolean hasCurveGeometry(int bLookForNonLinear = FALSE) const CPL_OVERRIDE;
186       virtual OGRErr addGeometry( const OGRGeometry * );
187       OGRErr addGeometryDirectly(OGRGeometry *poNewGeom);
188       int getNumGeometries() const;
189       OGRGeometry* getGeometryRef(int i);
190       const OGRGeometry* getGeometryRef(int i) const;
191
192       virtual OGRBoolean  IsEmpty() const CPL_OVERRIDE;
193       virtual void setCoordinateDimension( int nDimension ) CPL_OVERRIDE;
194       virtual void set3D( OGRBoolean bIs3D ) CPL_OVERRIDE;
195       virtual void setMeasured( OGRBoolean bIsMeasured ) CPL_OVERRIDE;
196       virtual void swapXY() CPL_OVERRIDE;
197       OGRErr removeGeometry( int iIndex, int bDelete = TRUE );
198   };
199
200-  The Triangulated Surface API is similar to Polyhedral Surface, and
201   the MultiPolygon class was tweaked slightly to include methods to run
202   which consisted of subgeometries of the form Triangle. (A
203   MultiPolygon is strictly a collection of Polygons). These methods are
204   internal to OGRMultiPolygon and cannot be accessed by a public user.
205   For instance, the ``OGRMultiPolygon::addGeometryDirectly`` method has
206   a check that the subgeometry added to it should be of the type
207   POLYGON. Rather than mess around with the existing function, a new
208   function has been written which does not implement this check -
209
210::
211
212   /************************************************************************/
213   /*                         _addGeometryDirectly()                       */
214   /*      Only to be used in conjunction with OGRTriangulatedSurface.     */
215   /*                        DO NOT USE IT ELSEWHERE.                      */
216   /************************************************************************/
217
218   OGRErr OGRMultiPolygon::_addGeometryDirectly( OGRGeometry * poNewGeom )
219   {
220       if ( wkbFlatten(poNewGeom->getGeometryType()) != wkbTriangle)
221           return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
222
223       if( poNewGeom->Is3D() && !Is3D() )
224           set3D(TRUE);
225
226       if( poNewGeom->IsMeasured() && !IsMeasured() )
227           setMeasured(TRUE);
228
229       if( !poNewGeom->Is3D() && Is3D() )
230           poNewGeom->set3D(TRUE);
231
232       if( !poNewGeom->IsMeasured() && IsMeasured() )
233           poNewGeom->setMeasured(TRUE);
234
235       OGRGeometry** papoNewGeoms = (OGRGeometry **) VSI_REALLOC_VERBOSE( papoGeoms,
236                                                sizeof(void*) * (nGeomCount+1) );
237       if( papoNewGeoms == NULL )
238           return OGRERR_FAILURE;
239
240       papoGeoms = papoNewGeoms;
241       papoGeoms[nGeomCount] = poNewGeom;
242       nGeomCount++;
243
244       return OGRERR_NONE;
245   }
246
247-  The Triangulated Surface API is as follows -
248
249::
250
251   class CPL_DLL OGRTriangulatedSurface : public OGRPolyhedralSurface
252   {
253     protected:
254   //! @cond Doxygen_Suppress
255       virtual OGRBoolean         isCompatibleSubType( OGRwkbGeometryType ) const CPL_OVERRIDE;
256       virtual const char*        getSubGeometryName() const CPL_OVERRIDE;
257       virtual OGRwkbGeometryType getSubGeometryType() const CPL_OVERRIDE;
258
259       virtual OGRPolyhedralSurfaceCastToMultiPolygon GetCasterToMultiPolygon() const CPL_OVERRIDE;
260       static OGRMultiPolygon* CastToMultiPolygonImpl(OGRPolyhedralSurface* poPS);
261   //! @endcond
262
263     public:
264       OGRTriangulatedSurface();
265       OGRTriangulatedSurface(const OGRTriangulatedSurface &other);
266       ~OGRTriangulatedSurface();
267
268       OGRTriangulatedSurface& operator=(const OGRTriangulatedSurface& other);
269       virtual const char *getGeometryName() const CPL_OVERRIDE;
270       virtual OGRwkbGeometryType getGeometryType() const CPL_OVERRIDE;
271
272       // IWks Interface
273       virtual OGRErr addGeometry( const OGRGeometry * ) CPL_OVERRIDE;
274
275       static OGRPolyhedralSurface* CastToPolyhedralSurface(OGRTriangulatedSurface* poTS);
276   };
277
278Geometry types
279--------------
280
281The new geometry WKB values can be seen as below -
282
283================= ==== ==== ==== ====
284Geometry Type     2D   Z    M    ZM
285================= ==== ==== ==== ====
286PolyhedralSurface 0015 1015 2015 3015
287TIN               0016 1016 2016 3016
288Triangle          0017 1017 2017 3017
289================= ==== ==== ==== ====
290
291Geometry conversions
292--------------------
293
294The OGRGeometryFactory::forceTo() and forceToMultiPolygon() methods have
295been enhanced to support conversions between the new geometry types, and
296towards multipolygon. Note that converting a TIN or a PolyhedralSurface
297into a MultiPolygon is semantically incorrect since a MultiPolygon is
298suppose to contain geometries in the same plane, but it might help when
299converting those new geometry types into a format that doesn't support
300them (and such conversion was for example implicitly done in the reading
301side of the shapefile driver previously)
302
303Changes in drivers
304------------------
305
306PostGIS
307~~~~~~~
308
309No changes done to the driver explicitly, but it has been ensured that
310PG <-> OGR compatibility has been maintained. PostGIS 3D functions work
311on OGR, simple scripts work, for example from
312``autotest/ogr/ogr_pg.py``, we have -
313
314::
315
316   wkt_list = ['POLYHEDRALSURFACE (((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((1 1 0,1 1 1,1 0 1,1 0 0,1 1 0)),((0 1 0,0 1 1,1 1 1,1 1 0,0 1 0)),((0 0 1,1 0 1,1 1 1,0 1 1,0 0 1)))',
317                   'TIN (((0 0 0,0 0 1,0 1 0,0 0 0)),((0 0 0,0 1 0,1 1 0,0 0 0)))',
318                   'TRIANGLE ((48 36 84,32 54 64,86 11 54,48 36 84))' ]
319
320   for i in range(0,3):
321           gdaltest.pg_ds.ExecuteSQL( "INSERT INTO zgeoms (field_no, wkb_geometry) VALUES (%d,GeomFromEWKT('%s'))" % ( i, wkt_list[i] ) )
322
323ShapeFile
324~~~~~~~~~
325
326Shapefiles have the concept of a "MultiPatch" object. The MultiPatch can
327be made of several parts, which are
328
329-  A TriangleStrip is a linked strip of triangles, where every vertex
330   (after the first two) completes a new triangle. A new triangle is
331   always formed by connecting the new vertex with its two immediate
332   predecessors.
333-  A TriangleFan is a linked fan of triangles, where every vertex (after
334   the first two) completes a new triangle. A new triangle is always
335   formed by connecting the new vertex with its immediate predecessor
336   and the first vertex of the part.
337-  Rings (outer ring, inner ring, first ring, "non-typed" ring) Up to
338   now multipatch were read as MultiPolygon. Now, in general, a
339   GeometryCollection will be returned, with zero or several TIN
340   corresponding to the TriangleStrip/TriangleFan and zero or one
341   MultiPolygon with all the rings. If there's only one TIN or one
342   MultiPolygon, it will be returned as a top-level geometry. The layer
343   type will be Unknown On writing, the SHPT layer creation option is
344   extended to recognize the MULTIPATCH value, and the current logic to
345   guess the shape type from the layer geometry type or the geometry
346   type of the first feature is extended to support MULTIPATCH. On a
347   MULTIPATCH layer, geometries of type TIN, POLYHEDRALSURFACE,
348   MULTIPOLYGON or GEOMETRYCOLLECTION (whose subgeometries are on of the
349   3 previous types) are accepted and converted to a MultiPatch object,
350   trying to use TriangleStrip and TriangleFan if the triangles are in
351   the expected order.
352
353FileGDB, OpenFileGDB
354~~~~~~~~~~~~~~~~~~~~
355
356The FileGDB format support the MultiPatch object as well, with one
357extension. There is a new type of part, which is made of several
358triangles whose organization is not TriangleStrip or TriangleFan. Both
359drivers have been upgraded to work like the ShapeFile driver on the
360reading side. On the writing side, the FileGDB driver will automatically
361write a MultiPatch if the layer geometry type is TIN or
362PolyhedralSurface. The layer option that existed before
363CREATE_MULTIPATCH=YES can still be used to force writing as MultiPatch
364
365GML
366~~~
367
368The GML driver has been modified for both input and output -> Triangle,
369PolyhedralSurface and TriangulatedSurface are capable of being
370read/written from/to a GML document. Sample examples include -
371
372::
373
374   'TRIANGLE ((0 0,0 1,0 1,0 0))' is parsed to -
375   '<gml:Triangle>
376       <gml:exterior>
377           <gml:LinearRing>
378               <gml:posList>0 0 0 1 0 1 0 0</gml:posList>
379           </gml:LinearRing>
380       </gml:exterior>
381   </gml:Triangle>'
382
383   <gml:PolyhedralSurface>
384      <gml:polygonPatches>
385          <gml:PolygonPatch>
386              <gml:exterior>
387                  <gml:LinearRing>
388                      <gml:posList srsDimension="3">1 2 3 4 5 6 7 8 9 1 2 3</gml:posList>
389                  </gml:LinearRing>
390              </gml:exterior>
391          </gml:PolygonPatch>
392          <gml:PolygonPatch>
393              <gml:exterior>
394                  <gml:LinearRing>
395                      <gml:posList srsDimension="3">10 11 12 13 14 15 16 17 18 10 11 12</gml:posList>
396                  </gml:LinearRing>
397              </gml:exterior>
398              <gml:interior>
399                  <gml:LinearRing>
400                      <gml:posList srsDimension="3">19 20 21 22 23 24 25 26 27 19 20 21</gml:posList>
401                  </gml:LinearRing>
402              </gml:interior>
403          </gml:PolygonPatch>
404      </gml:polygonPatches>
405   </gml:PolyhedralSurface>"""
406
407   gets parsed to 'POLYHEDRALSURFACE Z (((1 2 3,4 5 6,7 8 9,1 2 3)),((10 11 12,13 14 15,16 17 18,10 11 12),(19 20 21,22 23 24,25 26 27,19 20 21)))'
408
409   Each PolygonPatch/Patch corresponds to one Polygon in a PolyhedralSurface.
410
411   Finally, 'POLYHEDRALSURFACE EMPTY' parses to
412   '<gml:PolyhedralSurface>
413       <gml:polygonPatches>
414       </gml:polygonPatches>
415   </gml:PolyhedralSurface>'
416
417Note that on the writing side those geometries are only generated for a
418GML 3 output.
419
420DXF
421~~~
422
423The changes in the DXF driver include converting a PolyFaceMesh (a
424subtype of PolyLine) to PolyhedralSurface. This is illustrated by a bug
425on the GDAL trac -
426`https://trac.osgeo.org/gdal/ticket/6246 <https://trac.osgeo.org/gdal/ticket/6246>`__.
427A PolyFace Mesh consists of points defined initially using specific
428codes, then these points are described as part of a polygon (a polygon
429can have four points at the maximum). Reading the PolyFace Mesh is
430supported in OGR as of now, but write support for it as well (though not
431implemented by me in this changeset) should be possible as well now.
432
433GeoPackage
434~~~~~~~~~~
435
436The GeoPackage specification supports [Multi]Point, [Multi]LineString,
437[Multi]Polygon and GeometryCollection in its core. Curve geometry types
438are mentioned as a registered extension. But Triangle,
439PolyhedralSurface or TIN are not mentioned at all. However the
440GeoPackage geometry blob format being based on ISO WKB, support for the
441new geometry types did not really require new code. Hence we have kepts
442this possibility of reading/writing the 3 new geometry types, but with a
443warning emitted that a non-standard extension will be used on the
444writing side.
445
446Other drivers
447~~~~~~~~~~~~~
448
449The CSV, VRT, PGDump, SQLite (but not Spatialite) drivers support the
450new geometry types. A couple of drivers have been modified, so as not to
451crash on the writing side when being provided with the new geometry
452types. Besides the previously mentioned drivers, the following drivers
453have been verified to not crash (but potentially error out, or skip
454unrecognized geometries): MySQL, OCI, KML, LIBKML, GeoJSON, MapInfo
455
456Documentation
457-------------
458
459Using standard Doxygen documentation procedure.
460
461Compatibility
462-------------
463
464Many applications will not be able to properly deal with the new
465geometry types that may now be returned by some drivers. In GDAL 2.1,
466the new types were introduced mentioning that they might be returned by
467GDAL in the future. Code should either skip the new geometries, deal
468with them properly or use the OGR_G_ForceTo() function to convert to a
469geometry type it supports.
470
471Testing
472-------
473
474Very few changes have been made so that the existing autotest suite
475still passes. New geometry classes and conversion methods has been added
476to ogr_geom.py and ogr_gml_geom.py. Updated drivers have received new
477tests also.
478
479Implementation
480--------------
481
482Done by Avyav Kumar Singh, under the Google Summer of Code 2016 program,
483and fine tuned / extended / integrated by Even Rouault.
484
485The proposed implementation lies in the "gsoc-triangle-ps-tin-rebased"
486branch of the
487`https://github.com/rouault/gdal2/tree/gsoc-triangle-ps-tin-rebased <https://github.com/rouault/gdal2/tree/gsoc-triangle-ps-tin-rebased>`__
488repository.
489
490Voting history
491--------------
492
493+1 from JukkaR, DanielM, HowardB and EvenR
494