GIS Domain Knowledge Oracle — vibeSpatial

Predicate Edge Cases

Point on polygon boundary:

• within(point, polygon) = True (boundary contact allowed)
• contains(polygon, point) = True
• touches(point, polygon) = True (boundary-only intersection)
• contains_properly(polygon, point) = False (point is on boundary, not interior)
• intersects(point, polygon) = True

Point on polygon vertex: Same as point-on-boundary — vertices are part of the boundary.

Line endpoint touching polygon boundary:

• touches(line, polygon) = True (if only endpoint touches)
• intersects(line, polygon) = True
• crosses(line, polygon) = False (no interior penetration)

Identical geometries:

• within(a, a) = True
• contains(a, a) = True
• overlaps(a, a) = False for points/lines (same dim but identical = not "partial")
• equals(a, a) = True

Empty geometries:

• intersects(empty, anything) = False
• disjoint(empty, anything) = True
• All other predicates with empty input = False

Null geometries:

• Any predicate with null input = None (propagate null)

Predicate Symmetry

For asymmetric predicates: within(a, b) == contains(b, a).

2. Geometry Validity and Conventions

Ring Orientation (Winding)

• Exterior rings: Counter-Clockwise (CCW) = positive signed area
• Interior rings (holes): Clockwise (CW) = negative signed area
• Shoelace formula: signed_area = 0.5 * sum(x[i]*y[i+1] - x[i+1]*y[i])

A positive shoelace result means CCW. This is the OGC/ISO convention used by Shapely, PostGIS, and vibeSpatial.

Ring Closure

Rings MUST be closed: first coordinate == last coordinate. A triangle requires 4 coordinate pairs: [(0,0), (1,0), (0,1), (0,0)].

Minimum valid ring: 4 coordinates (3 unique + closing).

Polygon Validity Rules (OGC)

1. Exterior ring is CCW, holes are CW
2. All rings are closed
3. All rings are simple (no self-intersection)
4. Holes are contained within the exterior ring
5. No two rings cross (holes may touch exterior at a single point)
6. The interior is connected (no cuts that split it)

Invalid Geometry Examples

LineString Validity

• Minimum 2 points
• is_closed: first point == last point (forms a LinearRing)
• Self-intersection is allowed (unlike rings)

Point Validity

• Always valid unless null or empty
• Empty point: valid geometry with no coordinates

3. Measurement Semantics

Area

• Polygon: Absolute value of shoelace area (exterior - holes)
• MultiPolygon: Sum of individual polygon areas
• Point/LineString: 0.0 (dimensionally zero)
• Empty geometry: 0.0
• Null geometry: NaN (propagate)
• Units: Square units of the CRS (not geodetic!)

Length

• LineString: Sum of Euclidean segment lengths
• MultiLineString: Sum of part lengths
• Polygon: Perimeter (exterior ring + hole rings)
• Point: 0.0
• Empty: 0.0
• Null: NaN
• Units: Linear units of the CRS

Distance

• Euclidean: Minimum distance between closest points
• Same geometry: 0.0
• Disjoint: Positive value
• Intersecting: 0.0
• Null input: NaN
• Empty input: NaN (or infinity depending on convention)
• dwithin(a, b, d): distance(a, b) <= d

Centroid

• Polygon: Area-weighted centroid (not center of bounding box)
• LineString: Length-weighted midpoint
• Point: The point itself
• MultiPolygon: Area-weighted across parts
• Empty: Empty point
• Null: Null
• Centroid may be OUTSIDE the geometry (e.g., C-shaped polygon)

Bounds

• Returns (minx, miny, maxx, maxy)
• Empty geometry: (NaN, NaN, NaN, NaN)
• Null geometry: (NaN, NaN, NaN, NaN)
• Point: (x, y, x, y) (zero-extent box)

4. Constructive Operation Semantics

Intersection / Union / Difference / Symmetric Difference

These are set operations on the point sets of two geometries:

Return type depends on input and result dimension:

• Polygon intersect Polygon → Polygon, MultiPolygon, LineString, Point, or GeometryCollection
• Line intersect Line → Point, MultiPoint, LineString, or GeometryCollection
• Point intersect anything → Point or empty

Edge cases:

• intersection(a, a) = a (identity)
• union(a, a) = a
• difference(a, a) = empty
• intersection(a, disjoint_b) = empty
• Result may be mixed type (GeometryCollection) for non-trivial overlaps

Buffer

• buffer(geom, distance): All points within distance of geom
• buffer(geom, 0): For polygons, equivalent to make_valid (repair trick)
• buffer(point, r): Circle approximation (polygon with quad_segs segments per quarter)
• buffer(line, r): Corridor around the line
• buffer(polygon, -r): Erode inward (may produce empty if too thin)
• Negative buffer on point/line = empty
• cap_style: round (1), flat (2), square (3)
• join_style: round (1), mitre (2), bevel (3)

Clip (clip_by_rect)

• Fast path for axis-aligned rectangle clipping
• Sutherland-Hodgman algorithm
• Output preserves ring orientation
• Geometry fully inside clip box: returned unchanged
• Geometry fully outside clip box: empty geometry
• Partial overlap: clipped polygon (possibly multi-part)

Dissolve

• Group-by aggregation of geometries
• dissolve(by=column): Union all geometries sharing the same group key
• Interior boundaries between adjacent polygons are removed
• Produces MultiPolygon when unioned parts are disjoint

Make Valid

• Repairs invalid geometries following OGC rules
• Self-intersecting rings → split at intersection points → reassemble
• Inverted holes → reverse winding
• Unclosed rings → close by appending first coordinate
• Spike removal → collapse degenerate parts
• Output is always valid (may change geometry type)

5. Null and Empty Handling

This is the #1 source of subtle bugs. Null and empty are DIFFERENT.

Propagation Rules

• Metric ops (area, length, distance): Null → NaN, Empty → 0.0
• Predicates: Null → None, Empty → False (except disjoint)
• Constructive: Null → None, Empty → Empty (or context-dependent)
• Aggregation (dissolve): Nulls skipped, empties included

6. Coordinate System Considerations

vibeSpatial Assumptions

• All kernel computations are planar (Euclidean)
• CRS is stored at GeoSeries/GeoDataFrame level, not per-geometry
• Kernels are CRS-agnostic — they operate on raw (x, y) coordinates
• Distance, area, length are in CRS units (meters for UTM, degrees for WGS84)

When CRS Matters

• Area in degrees (WGS84/EPSG:4326): Meaningless for accurate measurement. Must project to equal-area CRS first. vibeSpatial does NOT do this automatically.
• Distance in degrees: 1 degree != constant meters. For accurate distance, project to local UTM or use geodetic formulas. vibeSpatial uses Euclidean distance only.
• Antimeridian (180/-180 longitude): Geometries crossing the antimeridian may have incorrectly computed bounds, area, and predicates. Split at antimeridian first.
• Polar regions: Euclidean assumptions break down near poles.

Coordinate Precision

• fp64 storage: ~15 significant decimal digits
• fp64 geographic: Sub-millimeter precision at any location on Earth
• fp32 geographic: ~1 meter precision (insufficient for parcel-level work)
• fp32 with centering: Sub-meter precision when coordinates are centered around the compute region

7. Degeneracy Corpus

vibeSpatial maintains a deterministic corpus of degenerate geometries for testing edge cases. When writing or verifying operations, test against these:

File: src/vibespatial/testing/degeneracy.py

8. Oracle Testing Contract

When verifying correctness, compare against Shapely (the reference implementation):

# Geometry equality
shapely_result.equals_exact(gpu_result, tolerance=1e-9)

# Float equality
math.isclose(gpu_value, shapely_value, rel_tol=1e-7, abs_tol=1e-9)

# Null handling
# Both null → equal
# One null, one not → NOT equal
# Both empty, same type → equal

Standard tolerances:

• rtol=1e-7 (relative) for fp64 compute
• atol=1e-9 (absolute) for fp64 compute
• For fp32 compute: rtol=1e-3 to 1e-4 depending on kernel class

9. Robustness Guarantees by Kernel Class

Exact predicates use adaptive precision:

1. Fast fp64 test with error bound check
2. If |orientation| <= error_bound: mark ambiguous
3. Re-evaluate ambiguous cases with higher precision (expansion arithmetic)
4. Error bound formula: error = (3.0 + 16.0 * eps) * eps * magnitude

10. "Does This Make Sense?" Quick Checks

Use these rules to sanity-check operation results:

• area(polygon) < 0: Wrong. Area is absolute value. If your kernel returns negative area, you have a winding bug.
• centroid outside polygon: Valid! C-shaped, U-shaped, and ring-shaped polygons can have exterior centroids.
• buffer(polygon, 0) != polygon: Valid if polygon was invalid. buffer(0) is a repair operation.
• intersection(a, b) larger than a or b: Wrong. Intersection is always <= both inputs.
• union(a, b) smaller than a or b: Wrong. Union is always >= both inputs.
• within(a, b) AND within(b, a): Means a equals b (geometrically).
• distance(a, b) = 0 but NOT intersects(a, b): Wrong. Zero distance implies intersection (or touching).
• touches(a, b) AND intersects(a, b): Always true. Touches implies intersects.
• touches(a, b) AND overlaps(a, b): Wrong. These are mutually exclusive.
• contains(a, b) AND disjoint(a, b): Wrong. Mutually exclusive.
• is_valid(make_valid(g)): Must always be True.
• area(union(a,b)) <= area(a) + area(b): Always true (equality when disjoint).