kenro 0.4.0

SpatiaLite-style spatial SQL for SQLite in pure Rust — PostGIS-compatible ST_ functions, GeoPackage R-tree, CRS transform, H3, MVT. Use via rusqlite, loadable extension, or WASM
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
# 3D geometry

> **Related:** [Function reference](functions.md) ·
> [Scope and semantics](scope.md) · [Quickstart](quickstart.md) ·
> [Transform accuracy](accuracy.md) · [WebAssembly hosts](wasm.md)

kenro computes in 2D: `geo_types` has two ordinates per coordinate and no room
for a third. And yet a height survives storage, indexing, every coordinate
transform, reprojection, and every derived geometry that can honestly keep one —
and where it cannot, the function says so rather than quietly handing back a
flattened result.

This page is how that works, in the order it happens to a geometry, with each
function table beside the semantics it needs. Everything asserted here was
**measured** against PostGIS 3.5 — the same image `scripts/golden/generate.sh`
generates the golden vectors from.

The short version:

| | |
|---|---|
| **read** | [pass-through](#3d-pass-through) — a Z stores, indexes, filters and reports; [surface collections](#surface-collections-polyhedralsurface-tin-triangle) — a POLYHEDRALSURFACE is read, measured and flattened, never decoded |
| **move** | [affine transforms](#3d-affine-transforms) and [reprojection](#st_transform) rewrite coordinates in the encoding, so a Z rides through |
| **derive** | [derived geometries](#derived-geometries-and-the-z) keep the heights of the vertices they reused, [interpolate](#interpolated-heights) between two of them — or refuse |
| **create** | [`ST_Force3D`](#creating-a-z) sets a height the caller names |
| **measure** | [the `ST_3D*` family](#3d-distance-and-predicates), which core PostGIS has without SFCGAL; [3D area and enclosed volume](#3d-area-and-enclosed-volume), which come from SFCGAL |

The one thing kenro will not do is **guess** a height: extrapolate one past the
end of a line, or average two that disagree.

---

## 3D pass-through

Decoding drops the Z and every encoder refuses a geometry that had one, rather
than silently writing 2D. These functions let a 3D column be *stored, indexed,
filtered and read* anyway — what a CityGML-style workflow needs even when the
analysis itself is planar.

The route in is the practical one: a GeoPackage written by GDAL or QGIS.
`ST_GeomFromGPB`, `ST_SetSRID` and `ST_AsGPB` carry the WKB payload across
byte-for-byte, so the Z survives — and so do surface collections, which those
three used to refuse because they validated by decoding. `ST_GeomFromWKB` still
re-encodes and so still flattens.

| Function | Returns | PostGIS | DuckDB Spatial | SpatiaLite | Notes |
|---|---|---|---|---|---|
| [`ST_HasZ(geom)`](https://postgis.net/docs/ST_HasZ.html) / [`ST_HasM(geom)`](https://postgis.net/docs/ST_HasM.html) | INTEGER | ✅ | ✅ | ✅ | Read from the encoding, not from kenro's (always 2D) decoded value |
| [`ST_NDims(geom)`](https://postgis.net/docs/ST_NDims.html) / [`ST_CoordDim(geom)`](https://postgis.net/docs/ST_CoordDim.html) | INTEGER | ✅ | ✅ | ✅ | 2, 3 or 4 — honestly. These used to answer a flat 2 |
| [`ST_Z(point)`](https://postgis.net/docs/ST_Z.html) / [`ST_M(point)`](https://postgis.net/docs/ST_M.html) | REAL / NULL | ✅ | ✅ | ✅ | NULL when the vertex has no such ordinate |
| [`ST_ZMin(geom)`](https://postgis.net/docs/ST_ZMin.html) / [`ST_ZMax(geom)`](https://postgis.net/docs/ST_ZMax.html) | REAL / NULL | ✅ | ❌ | ✅ | ⚠️ a 2D geometry answers **0**, not NULL — PostGIS derives these from a bbox whose Z slot is zero, and `WHERE ST_ZMax(g) > 100` should behave the same on both. NULL only for an empty geometry. Also accept `BOX3D(…)` text — see [`ST_3DExtent`](#st_3dextent) |

The **2D** functions stay planar on 3D input, which is the whole point of the
word: the R-tree columns, `ST_Distance`, `ST_Intersects`, `ST_Area` and the rest
answer about the geometry's shadow. Their 3D counterparts are a separate family,
[further down](#3d-distance-and-predicates).

What is still not here is a **3D geometry model** — there is no decoded value
with three ordinates anywhere in kenro, and every section below works on the
encoding instead. Two consequences follow, and they are the honest limits:

- **Most of the SFCGAL solid-modelling family is absent** — no
  `ST_3DIntersection`, `ST_3DUnion`, `ST_Extrude`, `ST_Tesselate`,
  `ST_MakeSolid` and no SOLID type. Those need CGAL's exact predicates and a
  topology model, which is a different project from a SQLite extension. The
  reference is no longer what blocks them: `postgis_sfcgal` is now loaded in
  the golden harness, and the two members that are *theorems* rather than
  library conventions have graduated — see
  [3D area and enclosed volume](#3d-area-and-enclosed-volume). See also
  [scope](scope.md#deliberately-out-of-scope).
- **A height is never guessed** — not extrapolated past the end of a line, not
  averaged between two that disagree. Interpolating *between two known heights*
  is fine and [is done](#interpolated-heights); inventing one is the line.

---

## Surface collections: POLYHEDRALSURFACE, TIN, TRIANGLE

kenro **reads** surface collections rather than decoding them. `geo_types`
has no variant for one, and a second geometry model would mean two
representations that can disagree, so these functions walk the encoding
instead: they answer structural questions, measure patch by patch, and hand
the whole thing to the 2D world through `ST_Force2D`. Two of them measure the
collection *as a solid shape* — [3D area and enclosed
volume](#3d-area-and-enclosed-volume) — and they too read the bytes.

The route in is a GeoPackage written by GDAL, QGIS or a CityGML importer —
the same route as [3D pass-through](#3d-pass-through). `ST_GeomFromWKB` will
not do it, because it re-encodes.

**Everything that needs a decoded geometry refuses surface input, loudly and
in one place.** A predicate or an overlay function raises with a message naming
`ST_Force2D`; silently flattening a building into overlapping faces would be
the same class of mistake as writing 2D where 3D went in. The exception is the
other function that reads the encoding rather than decoding it —
[`ST_Affine`](#3d-affine-transforms) — which transforms a surface collection
directly.

| Function | Returns | PostGIS | DuckDB Spatial | SpatiaLite | Notes |
|---|---|---|---|---|---|
| [`ST_NumPatches(geom)`](https://postgis.net/docs/ST_NumPatches.html) | INTEGER / NULL | ✅ | ❌ | ❌ | NULL for anything that is not a surface collection |
| [`ST_PatchN(geom, n)`](https://postgis.net/docs/ST_PatchN.html) | geometry / NULL | ✅ | ❌ | ❌ | Patch `n` as a 2D POLYGON, **1-based** like `ST_GeometryN`; Z readable via `ST_ZMin`/`ST_ZMax` |
| [`ST_GeometryType`](https://postgis.net/docs/ST_GeometryType.html) | TEXT | ✅ | ❌ | ❌ | `ST_PolyhedralSurface` / `ST_Tin` / `ST_Triangle` |
| [`ST_Dimension`](https://postgis.net/docs/ST_Dimension.html) | INTEGER | ✅ | ❌ | ❌ | 2, as PostGIS reports |
| [`ST_Area`](https://postgis.net/docs/ST_Area.html), [`ST_Perimeter`](https://postgis.net/docs/ST_Perimeter.html) | REAL | ✅ | ❌ | ❌ | Summed patch by patch — the **planar** sum PostGIS reports, not a 3D surface area (measured) |
| [`ST_NumGeometries`](https://postgis.net/docs/ST_NumGeometries.html), [`ST_IsEmpty`](https://postgis.net/docs/ST_IsEmpty.html) | INTEGER | ✅ | ❌ | ❌ | Patches count as members, so the R-tree triggers keep working |
| [`ST_MinX`](https://postgis.net/docs/ST_XMin.html) … [`ST_MaxY`](https://postgis.net/docs/ST_YMax.html), [`ST_ZMin`](https://postgis.net/docs/ST_ZMin.html), [`ST_ZMax`](https://postgis.net/docs/ST_ZMax.html), [`ST_HasZ`](https://postgis.net/docs/ST_HasZ.html), [`ST_NDims`](https://postgis.net/docs/ST_NDims.html) | | ✅ | ❌ | ❌ | Walked from the patches; a surface column stays indexable |
| [`ST_Force2D(geom)`](https://postgis.net/docs/ST_Force2D.html) | geometry | ✅ | ✅ | ❌ | → MULTIPOLYGON of the patches. ⚠️ a closed solid becomes **overlapping coplanar faces** — geometrically correct, visually surprising, and what PostGIS does |
| [`ST_IsClosed(geom)`](https://postgis.net/docs/ST_IsClosed.html) | INTEGER | ✅ | ❌ | ❌ | Is this a closed shell? Combinatorial, not geometric: every edge shared by exactly two patches, tested on the 3D coordinates |
| [`ST_3DArea(geom)`](https://postgis.net/docs/ST_3DArea.html) | REAL | ⚠️ SFCGAL | ❌ | ❌ | The **3D** surface area, unlike `ST_Area` above — [3D area and enclosed volume](#3d-area-and-enclosed-volume) |
| `kenro_volume(geom)` | REAL / NULL | ❌ | ❌ | ❌ | **kenro-only**, and the name is the point. See below |
| [`ST_Affine(…)`](https://postgis.net/docs/ST_Affine.html) | geometry | ✅ | ❌ | ❌ | Both arities transform a surface collection, patches and Z included, because they rewrite the encoding rather than decoding it — [3D affine transforms](#3d-affine-transforms) |
| `kenro_gpkg_extension_required(geom)` | TEXT / NULL | ❌ | ❌ | ❌ | **kenro-only.** See below |

`ST_GeomFromGML` also reads CityGML's surface wrappers — `gml:Solid`,
`gml:CompositeSurface`, `gml:Surface`, `gml:TriangulatedSurface`,
`gml:Triangle` — flattening them to a MULTIPOLYGON on the way in.

### The GeoPackage obligation

GeoPackage Annex F.1 makes an extended geometry type legal **only if the file
declares it**: one row in `gpkg_extensions` per (table, column).

```sql
INSERT INTO gpkg_extensions (table_name, column_name, extension_name, definition, scope)
VALUES ('buildings', 'geom', 'gpkg_geom_POLYHEDRALSURFACE',
        'http://www.geopackage.org/spec120/#extension_geometry_types', 'read-write');
```

…and `gpkg_geometry_columns.geometry_type_name` carries `POLYHEDRALSURFACE`
rather than a core type name.

Nothing enforces it — measured: GDAL 3.11 reads a POLYHEDRALSURFACE column
with the row missing, with a wrong `definition` URL, or with the wrong
`gpkg_geom_*` name, without a word. It warns only when the `extension_name`
is one it does not implement at all. So the row is a declaration to other
readers, not a gate, and a writer that omits it produces a file that works
and is out of spec.

kenro does not write that row either. It registers functions; it does not
manage schemas — the same reason SpatiaLite's `InitSpatialMetadata` is out of scope,
and a function with a side effect could no longer be `SQLITE_DETERMINISTIC`
and `SQLITE_INNOCUOUS`, which the GeoPackage triggers depend on. What it does
instead is **name the obligation** so it is detectable rather than folklore:

```sql
SELECT kenro_gpkg_extension_required(geom) FROM buildings LIMIT 1;
-- 'gpkg_geom_POLYHEDRALSURFACE', or NULL when no extension is needed
```

The name is deliberately not `GPKG_*`: the spec reserves the `gpkg` author
prefix for OGC-adopted extension *names*, and although it says nothing about
SQL function names, a `GPKG_` function would read as one the standard defines.

---

## 3D affine transforms

kenro computes in 2D, but a coordinate transform does not need a geometry
model: it needs each coordinate, once. So `ST_Affine` does not go through the
2D value every other function decodes into — it rewrites the coordinates
**in the encoding**, which means Z survives, `POLYHEDRALSURFACE` transforms,
and *placing a CityGML building into the world* works without kenro becoming
a 3D engine.

Three properties, each measured against PostGIS 3.5:

- **Z rides through the 2D form.** `ST_Affine(POINT Z (1 2 3), 2,0,0,2, 10,20)`
  is `POINT(12 24 3)` — PostGIS leaves the Z alone, and so does kenro.
  (Earlier kenro versions raised an error here instead.)
- **The 3D form cannot raise dimensionality.** On 2D input, `z` is taken as 0
  for the `x'`/`y'` rows and the `z'` row is discarded:
  `ST_Affine(POINT(1 2), 1,2,3, 4,5,6, 7,8,9, 10,20,30)` is `POINT(15 34)`,
  not a 3D point. Producing 3D from 2D is
  [`ST_Force3D`](#creating-a-z)'s job, and asking for it explicitly is the
  difference: a transform must not change dimensionality behind your back.
- **M is never mistaken for Z.** ISO dimension code 2 is XYM: three
  ordinates, none of them a height.
  `ST_Affine(POINT M (1 2 99), …, zoff := 30)` is `POINTM(11 22 99)`.

The 3D matrix is the upper 3×4 of a 4×4, row-major:

```text
x' = a·x + b·y + c·z + xoff
y' = d·x + e·y + f·z + yoff
z' = g·x + h·y + i·z + zoff
```

### Which functions take this path

Every function whose whole job is to move coordinates, and no others:

| On the encoding-level path | 2D only, deliberately |
|---|---|
| `ST_Affine` (both arities) | `ST_SnapToGrid` |
| `ST_Translate`, `ST_Scale` | `ST_ReducePrecision` |
| `ST_Rotate`, `ST_RotateZ` | |
| `ST_TransScale` | |
| `ST_FlipCoordinates` | |
| `ST_ShiftLongitude`, `ST_WrapX` | |
| `ST_Transform` | |

The split is not about effort — it is what PostGIS does. Measured on 3.5, the
left column all preserve Z and (except `ST_WrapX`) all accept a surface
collection. The right column **is not coordinate-wise there**: PostGIS drops
the vertices that collapse onto each other, so
`ST_SnapToGrid(LINESTRING(0 0,0.1 0.1,1 1,1.1 1.1), 1)` is `LINESTRING(0 0,1 1)`
and a fully-collapsing polygon is `POLYGON EMPTY`. kenro only rounds. Rewriting
coordinates in place would have handed those two 3D support while deepening
that divergence, so they keep to the 2D value — and now **raise an error on 3D
input instead of silently dropping the Z**, which is what they used to do.

`ST_WrapX` is on the left column but refuses surface collections, matching
PostGIS's own "Wrapping of PolyhedralSurface geometries is unsupported".

---

## Derived geometries and the Z

A coordinate transform is not the only thing that can keep a height. Anything
whose output coordinates *came from* its input's — `ST_Reverse`,
`ST_ExteriorRing`, `ST_Simplify`, `ST_ConvexHull` — could carry the Z along too,
and PostGIS does. kenro now does as well, and the rule is decided **per call
from the data** rather than from a list of function names:

1. No input carried a Z → nothing changes, and nothing costs.
2. Every coordinate of the result was a vertex of some input → the result is
   written with those heights.
3. A coordinate lies **between** two input vertices → its height is blended
   from them, by the **2D** distance ratio. That covers densifying, cutting and
   smoothing — see [interpolated heights](#interpolated-heights).
4. Some coordinate was **invented** with no single honest source — a buffer
   arc, an extrapolation past the end, a crossing where two surfaces disagree →
   **error**, naming `ST_Force2D`.

The last case is the point. The alternative kenro used to take was to return a
2D geometry without saying so.
`UPDATE buildings SET geom = ST_AsGPB(ST_Simplify(ST_GeomFromGPB(geom), 0.1))`
would flatten a whole table in silence. It now either keeps the heights or
refuses.

Deciding from the data also gets cases a hand-written list would have got
wrong: `ST_ConvexHull`'s output vertices *are* input vertices, so its Z
survives, while `ST_Segmentize`'s midpoints cannot — and the same
`ST_Intersection` call refuses or succeeds depending on whether the operands
actually cross.

| | |
|---|---|
| **Keeps the Z** (output reuses input vertices) | `ST_StartPoint` `ST_EndPoint` `ST_PointN` `ST_GeometryN` `ST_ExteriorRing` `ST_InteriorRingN` `ST_Boundary` `ST_Multi` `ST_Normalize` `ST_Reverse` `ST_ForcePolygonCW` `ST_ForcePolygonCCW` `ST_ForceRHR` `ST_AddPoint` `ST_SetPoint` `ST_RemovePoint` `ST_MakePolygon` `ST_MakeLine` `ST_LineFromMultiPoint` `ST_Points` `ST_RemoveRepeatedPoints` `ST_Simplify` `ST_SimplifyVW` `ST_ConvexHull` `ST_ConcaveHull` `ST_DelaunayTriangles` `ST_TriangulatePolygon` `ST_LineMerge` `ST_UnaryUnion` `ST_Subdivide` `ST_Split` `ST_MakeValid` `ST_Intersection` `ST_Difference` `ST_SymDifference` `ST_Union` — *when* the operation happens not to invent a vertex |
| **Refuses on 3D input** (the Z would have to be invented) | the same overlay and interpolation functions when they do: `ST_Segmentize` `ST_ChaikinSmoothing` `ST_LineSubstring` `ST_LineInterpolatePoint` `ST_LineInterpolatePoints` `ST_LineExtend` `ST_GeometricMedian`, the `ST_Union` **aggregate**, and any crossing overlay |
| **2D on purpose** (PostGIS answers 2D too — measured) | `ST_Centroid` `ST_PointOnSurface` `ST_Envelope` `ST_OrientedEnvelope` `ST_MinimumBoundingCircle` `ST_ClosestPoint` `ST_ShortestLine` `ST_LongestLine` `ST_Buffer` `ST_ClipByBox2D` `ST_AsMVTGeom` `ST_Force2D` |
| **2D on purpose** (a bounding box, whose corners are not input vertices) | `ST_Expand` `ST_BoundingDiagonal` `ST_MakeBox2D` `ST_Extent` |

Two details worth knowing:

- **A bounding box never borrows a neighbour's Z.** `ST_BoundingDiagonal` of a
  polygon whose (10 10) corner sits at z = 3 ends at that x and y but needs the
  box's zmax, 4. Answering 3 would be confidently wrong, so every box-shaped
  result stays 2D. ⚠️ PostGIS returns 3D for `ST_Expand` and
  `ST_BoundingDiagonal`; kenro's answer is 2D there.
- **Two heights at one plan position are ambiguous.** A vertical wall gives the
  same (x, y) two Z values, and there is no way to pick. Those coordinates
  count as "no honest Z", so the function refuses rather than choosing.

`ST_Project` is the one exception that asserts a Z for a coordinate no input
occupied: sliding a point along the ground does not change its elevation, which
is what PostGIS does too.

### Interpolated heights

Consecutive coordinates in one coordinate run — a linestring's vertices, a
ring's — are a segment, and the encoding says so without any geometry model
being involved. So a coordinate that lies on a segment gets the two ends'
heights blended by how far along it sits. That is what unblocks the
linear-referencing family:

| Function | Example, measured on PostGIS 3.5 |
|---|---|
| `ST_Segmentize` | `ST_Segmentize(LINESTRING Z (0 0 0,10 0 100), 2.5)` → `0 0 0, 2.5 0 25, 5 0 50, 7.5 0 75, 10 0 100`. Polygons densify **per ring** |
| `ST_LineSubstring` | `(…, 0.25, 0.75)` → `2.5 0 25, 7.5 0 75` |
| `ST_LineInterpolatePoint` / `…Points` | `(…, 0.5)` → `POINT(5 0 50)` |
| `ST_ChaikinSmoothing` | one pass → `0 0 0, 7.5 0 75, 12.5 0 75, 20 0 0` — 0.25/0.75 blends |
| `ST_Split` | the cut vertex takes the input's interpolated height; the **blade's own Z is ignored**, as in PostGIS |

**The fraction is 2D.** `LINESTRING Z (0 0 0,10 0 10,20 0 30)` has 2D lengths
10 and 10 but 3D lengths 14.14 and 22.36, so `ST_LineInterpolatePoint(…, 0.5)`
is the middle vertex — where only a 2D fraction puts it. (PostGIS's 3D-aware `ST_3DLineInterpolatePoint` is a separate
function, and it *is* implemented — in the
[3D metric family](#3d-distance-and-predicates), where it takes its fraction by
3D length.)

What still refuses, each for its own reason:

- ⚠️ **An overlay crossing.** Two surfaces crossing in plan view have no shared
  height. GEOS returns the **average** of the two — measured: operands at z = 0
  and z = 1000 give a crossing at 500, and swapping them returns the identical
  result, so it is a mean rather than a preference. kenro declines, because that
  number describes neither surface; the error names `ST_Force2D`. This is a
  deliberate divergence, and the widest kind: an error where PostGIS returns a
  value, which a ported query hits immediately. It is taken anyway because the
  alternative is a number that looks like data and measures nothing, and because
  kenro's rule is that it deviates loudly and documentedly rather than returning
  a quietly different result — an error satisfies that rule, reproducing an
  average these docs would have to call meaningless does not. Revisit if a real
  query wants the mean.
- ⚠️ **`ST_LineExtend`.** Its new vertex is *past* the last one, so no segment
  contains it. PostGIS extrapolates the final gradient (`… ,10 0 100,15 0 150`);
  kenro refuses rather than assuming a height keeps climbing.
- **`ST_ChaikinSmoothing` with more than one iteration.** The second pass works
  on chords between the first pass's points, which lie on no original segment.
  One pass succeeds, two refuse. (On a straight input every pass stays on the
  original segment, so it keeps working — correct rather than lucky.)
- **`ST_GeometricMedian`.** PostGIS computes it in 3D; there is nothing to
  interpolate from.
- **Two heights with equal claim**, from a vertical wall or a self-touching
  ring. Ambiguity is a positive statement that the answer is unknown.

---

## Creating a Z

`ST_Force3D` is the one function that *adds* a height rather than carrying or
computing one, and it is a fair amount less exotic than it sounds: the XYZ
writer built to carry heights across a derived geometry already emits ISO XYZ
type codes, so all `ST_Force3D` supplies is a constant instead of a lookup. No
decoded 3D geometry model is involved.

| Function | Returns | PostGIS | DuckDB Spatial | SpatiaLite | Notes |
|---|---|---|---|---|---|
| [`ST_Force3D(geom [, zvalue])`](https://postgis.net/docs/ST_Force_3D.html) | geometry | ✅ | ❌ | ⚠️ named `CastToXYZ` | Every coordinate gains a Z, default 0. Works on every type, including collections |
| [`ST_Force3DZ(geom [, zvalue])`](https://postgis.net/docs/ST_Force_3DZ.html) | geometry | ✅ | ❌ | ❌ | PostGIS's alias for the same function |
| [`ST_MakePoint(x, y, z)`](https://postgis.net/docs/ST_MakePoint.html) | geometry | ✅ | ✅ | ⚠️ named `MakePointZ` | The four-argument XYZM form is not implemented — kenro cannot write an M |

Measured on PostGIS 3.5, and matched:

- `ST_Force3D(POINT(1 2))` is `POINT(1 2 0)`; with `zvalue = 7`, `POINT(1 2 7)`.
- **An existing Z is never overwritten.** `ST_Force3D(POINT Z (1 2 3), 7)` is
  `POINT(1 2 3)` — the argument fills gaps rather than setting heights. Use
  `ST_Affine`'s 3D form or `ST_Translate` to *change* a height.
- **XYM loses its M**: `ST_Force3D(POINT M (1 2 99))` is `POINT(1 2 0)`. The
  result is XYZ, and kenro has no XYZM writer.
- An empty geometry has no ordinates, so it comes back unchanged.

⚠️ A **surface collection with no Z in its type code** is refused rather than
raised: adding an ordinate there means rebuilding the nested patch encoding,
which this writer does not do. Surfaces that already carry a Z — which is all of
them in practice — pass straight through.

### `ST_Transform`

Reprojection is coordinate-wise, so it takes the encoding-level path as well —
which matters more here than anywhere else, because reprojecting is the
operation a 3D city model needs most and it used to refuse 3D outright.
Measured on PostGIS 3.5:

- `ST_Transform(POINT Z (139.7 35.7 100), 32654)` moves x and y and returns
  `z = 100` untouched.
- A `POLYHEDRALSURFACE Z` comes back a `POLYHEDRALSURFACE`, patch structure and
  roof heights intact, rather than being rejected.

The Z is not merely carried past the projection: proj4rs takes `(x, y, z)` per
coordinate, so a datum shift routed through geocentric coordinates *reads* the
height, exactly as PROJ does. Same-datum pairs (`4326 → 32654` among them)
leave it exactly as it was — verified by moving the height and watching x and y
not move.

The byte-level I/O functions came along for the ride, because a reprojected
building has to be storable: `ST_SetSRID`, `ST_GeomFromGPB`, `ST_AsGPB` and
`ST_SRID` all used to validate by *decoding*, which refused surface
collections — breaking the pass-through promise for exactly the geometries it
was written for. They now walk the encoding instead, and `ST_AsGPB` builds a
surface's R-tree envelope from its patches.

### `ST_3DExtent`

⚠️ PostGIS returns its `box3d` type. SQLite has no such type, and kenro cannot
write a 3D geometry to stand in for one, so `ST_3DExtent` returns **the text
PostGIS renders a box3d as**: `BOX3D(minx miny minz,maxx maxy maxz)`. The
digits are Rust's shortest round-trip, not PostGIS's — PostGIS renders a box3d
through the server's `extra_float_digits`, so its own output is not a fixed
string either.

**The six box accessors read that text back.** In PostGIS `ST_XMin`,
`ST_XMax`, `ST_YMin`, `ST_YMax`, `ST_ZMin` and `ST_ZMax` have exactly one
overload each and its argument type *is* `box3d` — a geometry only reaches
them through an implicit cast. SQLite has neither the type nor the cast, so
kenro's six take the union directly: a geometry BLOB, or the box text.

```sql
SELECT ST_XMin(b), ST_YMin(b), ST_ZMin(b),
       ST_XMax(b), ST_YMax(b), ST_ZMax(b)
FROM (SELECT ST_3DExtent(geom) AS b FROM buildings);
```

`BOX(minx miny,maxx maxy)` — the 2D spelling PostGIS's own `ST_Extent`
renders — parses too, with Z = 0, as do lower case and surrounding
whitespace. PostGIS's parser refuses all three (it is `sscanf`, so it is
case-sensitive), but in SQLite there is no cast to route a value through and
the typed string is the only way in. In the other direction kenro is
*stricter* than PostGIS about the same parser's accidents: a trailing junk
tail, a missing close paren and mismatched corner dimensions are errors here
and silently accepted there. `src/functions/box3d.rs` carries the measured
table, and `tests/golden/box_text.jsonl` pins every row of it.

For a whole-column box the SQLite aggregates over kenro's scalars remain the
direct route, and skip the round trip through text:

```sql
SELECT min(ST_MinX(geom)), min(ST_MinY(geom)), min(ST_ZMin(geom)),
       max(ST_MaxX(geom)), max(ST_MaxY(geom)), max(ST_ZMax(geom))
FROM buildings;
```

A 2D row contributes Z = 0 rather than nothing, following `ST_ZMin`/`ST_ZMax`
(PostGIS answers `BOX3D(0 0 0,5 5 0)` for `LINESTRING(0 0,5 5)`). Empty
geometries contribute nothing; a zero-row or all-empty group is NULL.

---

## 3D distance and predicates

The nine functions core PostGIS has **without** SFCGAL. What makes them possible
without a 3D geometry model is a single measurement: a POLYHEDRALSURFACE is
treated as a *set of faces*, never as a volume — a point at the dead centre of a
closed unit cube gives `ST_3DIntersects = false`. So there is no point-in-solid
test, no shell orientation and no topology, only primitives against primitives.

| Function | Returns | PostGIS | DuckDB Spatial | SpatiaLite | Notes |
|---|---|---|---|---|---|
| [`ST_3DDistance(a, b)`](https://postgis.net/docs/ST_3DDistance.html) | REAL / NULL | ✅ | ❌ | ✅ | NULL for an empty operand |
| [`ST_3DDWithin(a, b, d)`](https://postgis.net/docs/ST_3DDWithin.html) | INTEGER | ✅ | ❌ | ✅ | `false` for an empty operand, not NULL |
| [`ST_3DDFullyWithin(a, b, d)`](https://postgis.net/docs/ST_3DDFullyWithin.html) | INTEGER | ✅ | ❌ | ✅ | the maximum distance is at most `d` |
| [`ST_3DMaxDistance(a, b)`](https://postgis.net/docs/ST_3DMaxDistance.html) | REAL / NULL | ✅ | ❌ | ✅ | **vertex to vertex** — measured |
| [`ST_3DIntersects(a, b)`](https://postgis.net/docs/ST_3DIntersects.html) | INTEGER | ✅ | ❌ | ✅ | faces are filled, solids are not |
| [`ST_3DClosestPoint(a, b)`](https://postgis.net/docs/ST_3DClosestPoint.html) | geometry / NULL | ✅ | ❌ | ✅ | the point on `a` |
| [`ST_3DShortestLine(a, b)`](https://postgis.net/docs/ST_3DShortestLine.html) | geometry / NULL | ✅ | ❌ | ✅ | a 3D LINESTRING |
| [`ST_3DLongestLine(a, b)`](https://postgis.net/docs/ST_3DLongestLine.html) | geometry / NULL | ✅ | ❌ | ✅ | the vertex pair |
| [`ST_3DLineInterpolatePoint(line, f)`](https://postgis.net/docs/ST_3DLineInterpolatePoint.html) | geometry | ✅ | ❌ | ❌ | the fraction is by **3D** length, unlike the 2D sibling |

Points, linestrings, polygons, MULTI\* and the surface collections
(POLYHEDRALSURFACE, TIN, TRIANGLE) are all accepted — the CityGML shapes
included, which core PostGIS supports too. 196 golden vectors from the reference
cover every type pair; the suite is `tests/golden/threed.jsonl`.

Why they differ from their 2D namesakes, in one example: two lines that cross in
plan but sit 4 apart in height give `ST_Intersects = true` and
`ST_3DIntersects = false`.

Three behaviours are copied deliberately, because each is wrong by default:

- **A missing Z means "any value", not zero.** PostGIS says so in a notice, and
  `ST_3DDistance(POINT Z (0 0 10), POINT(0 0))` is **0** — the Z-less operand
  behaves as a vertical line. kenro delegates to the 2D functions there, which is
  exactly equivalent.
- **`ST_3DMaxDistance` measures vertices**, not interiors: a 10×10 face against
  its own corner is √200.
- **A non-planar ring is triangulated**, not flattened to a best-fit plane. A
  ring whose corners sit at z = 0, 0, 0, 10 answers 90.27735042633894 against a
  point at z = 100, where a planar ring answers exactly 100. kenro fans from the
  ring's first vertex; on that case every candidate diagonal gives the same
  minimum.

⚠️ **Two documented divergences, both where the reference contradicts itself.**

- **An empty operand.** `ST_3DShortestLine(POINT Z (0 0 0), LINESTRING EMPTY)`
  returns `LINESTRING(0 0 0,0 4.63557111106545e-310 0)` in PostGIS 3.5 — that
  subnormal is uninitialised memory. kenro returns NULL, which is what its own
  `ST_3DDistance` answers for the same pair.
- **A point at the exact centre of a coplanar face.**
  `ST_3DDistance(POINT Z (5 5 0), POLYGON Z ((0 0 0,10 0 0,10 10 0,0 10 0,…)))`
  is 7.0710678118654755 there — the distance to a corner — while
  `ST_3DClosestPoint` on the same pair returns `POINT(5 5 0)`, i.e. distance 0,
  and `ST_3DIntersects` returns `false`. Three answers that cannot all be right.
  Nearby interior coplanar points all answer 0 and `true`. kenro answers 0 and
  `true`, which is what the reference's own `ST_3DClosestPoint` implies.

---

## 3D area and enclosed volume

Two functions from SFCGAL's column, and the only two whose answer is a
**theorem** rather than a library's convention: the area of a surface is a sum
of cross products, and the volume a closed shell encloses is the divergence
theorem. Neither needs an exact-arithmetic kernel, a topology model or a new
dependency — which is why they are here while `ST_3DIntersection` and the rest
of that family are not.

The reference for both is `postgis_sfcgal` 1.3.8 in the same image every other
suite uses; `scripts/golden/generate.sh` loads it, and the vectors are
`tests/golden/threed_sfcgal.jsonl`.

| Function | Returns | PostGIS | DuckDB Spatial | SpatiaLite | Notes |
|---|---|---|---|---|---|
| [`ST_3DArea(geom)`](https://postgis.net/docs/ST_3DArea.html) | REAL | ⚠️ SFCGAL | ❌ | ❌ | The area of the faces measured in 3D. 0 for non-areal or empty, never NULL |
| `kenro_volume(geom)` | REAL / NULL | ❌ | ❌ | ❌ | **kenro-only.** The **signed** volume a closed, consistently oriented shell encloses. NULL for anything that is not a surface collection; a loud error for an open or inconsistent one |

### Why one of them is not called `ST_Volume`

The measurement that decided it, on a closed box of 3.3 × 1.7 × 3.6:

| input | `ST_3DArea` | `ST_Volume` |
|---|---|---|
| the box as POLYHEDRALSURFACE | 47.22 | **0** |
| the box as TIN | 47.22 | **0** |
| `ST_MakeSolid` of the POLYHEDRALSURFACE | **0** | 20.196 |

A surface encloses nothing, so SFCGAL's `ST_Volume` answers **0** for a closed
POLYHEDRALSURFACE; the volume appears only once the same coordinates are wrapped
as a SOLID, and the two functions are near-complements across the split.
[kenro has no SOLID type](scope.md#deliberately-out-of-scope), so a kenro
`ST_Volume` returning 20.196 for the shell would be a silently different result
under a shared name — the one thing kenro's naming rule forbids. `kenro_volume`
wears its own name instead, exactly as
[`kenro_gpkg_extension_required`](#the-geopackage-obligation) does.

`ST_3DArea` has no such problem: kenro's answer *is* SFCGAL's answer for every
encoding kenro has, so it wears the PostGIS name.

### What `kenro_volume` refuses, and why

`Σ (1/6) a·(b × c)` over the triangulated faces means something only when the
faces bound a region and all agree which side is out. Both conditions are one
test: **in a closed, consistently oriented shell every directed edge occurs
exactly once, and its reverse exactly once.**

- an edge with no partner → an **open shell** → `kenro: kenro_volume: not a
  closed shell …`
- a directed edge seen twice → a **flipped face** → `kenro: kenro_volume: the
  patches are not consistently oriented …`

SFCGAL refuses both too — `ST_MakeSolid` on an open box gives "shell 0 is not
closed", and one reversed ring makes it refuse before measuring anything at all,
`ST_3DArea` included.

That gate is **stricter than `ST_IsClosed`**, which counts edges *undirected*: a
box with one face reversed is `ST_IsClosed = true` and is still refused here.

**The sign is kept.** SFCGAL's solid volume is signed by shell orientation —
reversing every face turns 20.196 into −20.196, and an irregular tetrahedron's
−8.9155 into +8.9155 — so `kenro_volume` reproduces the sign rather than
throwing it away. Wrap in `abs()` for a magnitude; keep it to tell an
outward-facing shell from an inward-facing one.

```sql
SELECT abs(kenro_volume(geom)), ST_3DArea(geom) FROM buildings WHERE ST_IsClosed(geom);
```

### Where kenro is less strict than the reference

SFCGAL *validates* a surface before measuring it. kenro measures what it was
given, so two inputs answer where SFCGAL raises. Both are recorded as
divergence vectors, and neither is a case where both produce a number and the
numbers differ:

- **A flipped ring.** `ST_3DArea` answers 47.22 for the box with one face
  reversed, because a face's area does not depend on which way its ring runs.
  (`kenro_volume` does refuse it — there the orientation is the whole question.)
- **A degenerate or non-planar ring.** Three collinear vertices give 0, the area
  such a ring actually has, where SFCGAL says "Polygon is invalid"; a non-planar
  ring gets the area of its best-fit planar projection, where SFCGAL says "points
  don't lie in the same plane". kenro has no 3D validity checker, and inventing
  one to reproduce a refusal would be a larger claim than the measurement
  supports.

Holes subtract, and non-convex rings are exact: an L-shaped hexagon answers 3,
which a fan of triangle magnitudes would overcount. Rings arrive from the
encoding with their ordinal, so a polygon's vector areas are summed before the
magnitude is taken — an interior ring is wound against its shell, so the
subtraction *is* the sum.

---