geometry-simple-0.1.1.0: docs/GEOS-DIFFERENCES.md

# Simple Features and GEOS

The target is the planar Simple Features core for seven geometry families.
GEOS is an independent comparison reference. It does not define every Haskell
API or output convention in this package.

## Deliberate differences

| Behavior | geometry-simple | GEOS 3.13.1 |
| --- | --- | --- |
| Constructed planar coordinates | Hulls, centroids, representative points, overlays, and buffers use XY, including empty results. | Some operations retain, interpolate, or discard Z/M according to operation-specific rules. |
| Polygon construction | Exterior rings run counterclockwise; holes run clockwise. | Hulls and many constructed polygon exteriors run clockwise. |
| Hull ordering | The first vertex and line endpoints follow lexicographic XY order. | Polygon starts use Y then X; two-point hulls can retain input order. |
| Overlay line components | Join consecutive edges through vertices with two neighbors. | Can retain separate lines at source vertices and intersection nodes. |
| Representative points | Use the first available polygon, interior line vertex, endpoint, or point, in that order. Empty components are skipped. | Selection can depend on centroid distance, interval width, and empty members. |
| Point observers | Preserve the stored layout and every ordinate, including NaN Z/M. | Extracted points can lose dimensions whose ordinate is NaN. |
| Mixed WKT collections | Each child carries its own dimension tag; the parent has none. | A parent tag can conflict with a child and make the writer's output unreadable. |
| Bare MULTIPOINT with EMPTY | Accepts empty members before, between, or after bare coordinates, as emitted by DuckDB. Nonempty members must use one spelling throughout. | Requires parenthesized nonempty members when EMPTY is present. |
| WKT numbers | Shortest scientific notation that decodes to the same `Double`. | Decimal formatting differs. |
| Numerical reductions | Compensated centroid sums retain small contributions during cancellation. | Evaluation order and final floating-point digits can differ. |

OGC 06-103r4 section 6.1.2.5 distinguishes observers, which retain stored Z/M,
from map operations, which ignore them when computing new geometry.
Section 6.1.11.1 describes exterior rings as counterclockwise when viewed from
the top, with opposite winding for holes. We use that convention in the XY
plane. Lexicographic start vertices are our deterministic convention; the
standard does not prescribe them.

The SQL test annex also permits either polygon winding when checking its
expected answers. A different winding alone does not establish a topology bug.

Section 6.1.2.4 defines overlays by their point sets. It does not require a
particular grouping into line components. Sections 6.1.10.2 and 6.1.13.2 require
`PointOnSurface` to return a point on the surface. They do not prescribe which
point to select. The Haskell method also supports points, curves, and mixed
collections, preferring a nonempty component of the highest dimension.

Sources: [OGC common architecture](https://docs.ogc.org/is/06-103r4/06-103r4.pdf),
[OGC SQL test examples](https://docs.ogc.org/is/06-104r4/06-104r4.pdf).

## Geometry contracts

- Point, line, and polygon empties retain their family's topological dimension.
  A collection with no atomic members has dimension -1.
- Coordinate dimensions combine stored member metadata. XYZ and XYM members
  together have coordinate dimension 3, while both Z and M flags are present.
- Empty multi-geometries and empty collections have no stored coordinate layout.
  Writers inherit the containing collection's output tag, or use XY at the top
  level. An empty container does not force a collection to use mixed layouts.
- Polygon area subtracts holes regardless of input winding. Area and perimeter
  close rings supplied directly through Haskell constructors.
- `geometryLength` includes lines and polygon boundaries. `curveLength` counts
  only lines; `perimeter` counts only polygon rings.
- A centroid weights polygons by area, then falls back to segment lengths, then
  points. Each collapsed line or ring contributes its first coordinate once.
  Lower-dimensional components do not affect a higher-dimensional centroid.
- Envelopes use the bounding rectangle specified in section 6.1.2.2. Empty
  input gives an empty point; one XY location gives a point. Flat bounds retain
  repeated rectangle corners and can form a topologically invalid polygon.
- Line boundaries use the mod-2 endpoint rule. `boundary` returns `Nothing` for
  a geometry collection. Stored polygon rings retain their coordinate layouts.
- `relate` returns a nine-character DE-9IM matrix. `relatePattern` accepts
  `T`, `F`, `*`, `0`, `1`, and `2`; invalid patterns return `False`.
  Distance to an empty geometry is NaN.
- Buffers have round caps and joins, with eight segments per quadrant by
  default. Negative distances erode polygons. For points and lines, they give
  an empty polygon.
  Zero distance extracts polygonal regions. Invalid input can lose regions,
  such as one lobe of a self-crossing bowtie. Circular arcs are approximations.
- Overlay line results join through vertices with exactly two neighbors.
  They stop at endpoints and branches. The point set is preserved, but line
  component counts, component order, and polygon hole order can differ from GEOS.
- `pointOnSurface` returns an XY point on a nonempty component. Polygons use
  a horizontal interior interval. Lines prefer stored interior vertices over
  endpoints. A collection of empty components gives an empty XY point.
- Measured queries select points and curve portions using M. Polygon queries
  select their boundary positions, as permitted by the implementation-defined
  surface rule. Empty input gives `Nothing`; no match gives an empty point.

## Numerical limits

Planar operations require finite XY coordinates. Measurements use `Double`.
Polygon cross products use coordinates relative to a ring vertex. Centroid
sums keep polygon positions separate from local moments and compensate for
rounding during addition. Products can still overflow or lose precision.

Centroids can overflow for polygon spans above about 1e100 units or line spans
above about 1e150 units. Underflow can occur below about 1e-100 and 1e-150,
respectively. Compensation cannot recover rounding already lost in products.

Hull orientation tests use a bounded `Double` calculation with an exact
`Rational` fallback. Topology uses exact rational intersections and rounds
constructed output coordinates to `Double`. Exact bounding-box indexes prune
segment pairs, validity checks, point locations, and distance candidates.
Representative points use exact scanline crossings. Their rounded X ordinate
must stay within the selected interval. A shell vertex is used when no interval
contains a representable point. Distance retains exact scaling until the final
conversion, including subnormal projected distances.

Buffers use scaled norms when intermediate products overflow or underflow.
An offset outside the finite `Double` range returns `Left CoordinateOverflow`.
These rules can differ from GEOS at extreme coordinate scales.

Repeated winding queries use aggregated crossing counts. Arrangements in which
many segment bounds overlap can still take quadratic time, even when the shapes
do not intersect, such as slanted interleaved combs. Predicates can reject incompatible dimensions
or bounds before constructing a full relation matrix.

### Overlay precision

Overlays and buffers first compute an exact arrangement and round its result
to `Double`. Buffer offsets approximate circular arcs with floating-point coordinates.
Rounding can move all vertices of a small ring onto one line. The operation
removes such a ring and does not retry because of it.
A valid result keeps its coordinates. If rounding makes it invalid, the operation
retries with vertex and segment snapping. Separate groups of overlapping input
bounds use separate tolerances. A distant component therefore does not set the
precision of a local operation. Positive buffers expand these bounds by their
distance before grouping, because offsets from disjoint inputs can overlap.

The first tolerance is the group's largest absolute ordinate divided by 10^12,
with a floor of the smallest positive `Double`. Five attempts increase this
tolerance tenfold, each starting from the original inputs. Nearby vertices and
intersection nodes share coordinates. Edges also snap to nearby vertices.
Narrow regions can collapse. These functions return
`Either TopologyException Geometry`. If all precision attempts fail, they return
`Left PrecisionFailure`. Offsets beyond the finite `Double` range return
`Left CoordinateOverflow`. `Left OpenBoundary` and `Left UncontainedHole` indicate
a library defect. `Data.Geometry.SimpleFeatures` exports these
constructors. Applications can handle failures without catching exceptions.
The combined output is checked again after separate groups are processed.

The retry count and tolerance schedule follow
[GEOS OverlayNGRobust](https://github.com/libgeos/geos/blob/3.13.1/include/geos/operation/overlayng/OverlayNGRobust.h).
This implementation retains exact noding and deterministic XY representatives.
Buffers use this same retry schedule on their source coordinates; this differs
from GEOS buffer-specific precision reduction.
It does not reproduce GEOS's additional self-union and precision-grid attempts.
It can therefore report a precision failure before GEOS exhausts its retry
strategies. The two implementations do not have the same failure behavior.
Valid results can retain tiny and thin regions that GEOS discards. At subnormal scales,
GEOS's floating-point validity checks can also disagree with exact orientation.
The precision tests include an independent rational check for such a region.

## Format rules

Untagged WKT infers XY, XYZ, or XYZM from two, three, or four ordinates. XYM
requires M. Untagged collection children infer their layouts independently.
In multi-geometries, empties before the first coordinate remain XY; later
empties use the inferred layout. Explicit parent tags require matching child
layouts. Writers tag collections when all members share one output layout,
including nested collections. Empty containers inherit the containing layout;
standalone empty containers use XY.

Mixed-layout collection WKT omits the parent tag and retains each child's tag.
This extends the OGC grammar in section 7. GEOS accepts this form. DuckDB's
WKT reader and `ST_AsText` require one layout across all members. DuckDB can
store mixed-layout WKB through `ST_GeomFromWKB` and return it through
`ST_AsWKB`, but `ST_AsText` rejects that geometry. Use a common layout when
DuckDB interchange includes WKT.

The WKT decoder accepts attached tags such as `POINTZ`, both multipoint
syntaxes, signed numbers, fractions, and exponents. It also accepts `EMPTY`
with bare MULTIPOINT coordinates, including DuckDB's `ST_AsText` output. This
syntax is broader than GEOS 3.13.1. Mixing parenthesized and bare nonempty
members remains an error. Whitespace between ordinates
is required: space, tab, CR, or LF. The decoder rounds each decimal ordinate
to the nearest `Double`. Underflow produces signed zero; overflow produces infinity.

Both codecs accept nonfinite ordinates for storage. NaN in both WKB point XY
ordinates denotes an empty point. Nonfinite values do not have the finite-bit
round-trip guarantee. Use `EmptyPoint` for an empty point in Haskell.

Lines must have zero or at least two coordinates. Rings must have zero or at
least three coordinates and close in XY. An empty exterior cannot contain a
nonempty hole. Writers normalize polygons with only empty rings to an empty
polygon, so empty-hole counts can change after writing. These construction
checks do not establish valid topology.

WKB rejects unknown type tags, incorrect child families, and counts exceeding
the remaining input before allocation. Counts must fit in 32 bits. The codecs
have no nesting limit, and WKT has no count limit. EWKB, EWKT, and embedded SRIDs
are outside the package's scope.