fix: webmercator panic (shifting quadrants) (#35)

* add test for deviation in webmercatorquad
* fix shifting quadrants in pointindex because of int division rounding
* add DeviationStats and IsQuadTree

PDOK-16504

intgeom/extent.go

@@ -6,20 +6,20 @@ import (
 
 // Extenter represents an interface that returns a boundbox.
 type Extenter interface {
-	Extent() (extent [4]int64)
+	Extent() (extent [4]M)
 }
 
 // MinMaxer is a wrapper for an Extent that gets min/max of the extent
 type MinMaxer interface {
-	MinX() int64
-	MinY() int64
-	MaxX() int64
-	MaxY() int64
+	MinX() M
+	MinY() M
+	MaxX() M
+	MaxY() M
 }
 
 // Extent represents the minx, miny, maxx and maxy
 // A nil extent represents the whole universe.
-type Extent [4]int64
+type Extent [4]M
 
 func (e Extent) ToGeomExtent() geom.Extent {
 	return geom.Extent{
@@ -43,8 +43,8 @@ func FromGeomExtent(e geom.Extent) Extent {
 
 // Vertices return the vertices of the Bounding Box. The vertices are ordered in the following manner.
 // (minx,miny), (maxx,miny), (maxx,maxy), (minx,maxy)
-func (e Extent) Vertices() [][2]int64 {
-	return [][2]int64{
+func (e Extent) Vertices() [][2]M {
+	return [][2]M{
 		{e.MinX(), e.MinY()},
 		{e.MaxX(), e.MinY()},
 		{e.MaxX(), e.MaxY()},
@@ -54,15 +54,15 @@ func (e Extent) Vertices() [][2]int64 {
 
 // ClockwiseFunc returns whether the set of points should be considered clockwise or counterclockwise.
 // The last point is not the same as the first point, and the function should connect these points as needed.
-type ClockwiseFunc func(...[2]int64) bool
+type ClockwiseFunc func(...[2]M) bool
 
 // Edges returns the clockwise order of the edges that make up the extent.
-func (e Extent) Edges(cwfn ClockwiseFunc) [][2][2]int64 {
+func (e Extent) Edges(cwfn ClockwiseFunc) [][2][2]M {
 	v := e.Vertices()
 	if cwfn != nil && !cwfn(v...) {
 		v[0], v[1], v[2], v[3] = v[3], v[2], v[1], v[0]
 	}
-	return [][2][2]int64{
+	return [][2][2]M{
 		{v[0], v[1]},
 		{v[1], v[2]},
 		{v[2], v[3]},
@@ -71,36 +71,36 @@ func (e Extent) Edges(cwfn ClockwiseFunc) [][2][2]int64 {
 }
 
 // MaxX is the larger of the x values.
-func (e Extent) MaxX() int64 {
+func (e Extent) MaxX() M {
 	return e[2]
 }
 
 // MinX  is the smaller of the x values.
-func (e Extent) MinX() int64 {
+func (e Extent) MinX() M {
 	return e[0]
 }
 
 // MaxY is the larger of the y values.
-func (e Extent) MaxY() int64 {
+func (e Extent) MaxY() M {
 	return e[3]
 }
 
 // MinY is the smaller of the y values.
-func (e Extent) MinY() int64 {
+func (e Extent) MinY() M {
 	return e[1]
 }
 
 // XSpan is the distance of the Extent in X
-func (e Extent) XSpan() int64 {
+func (e Extent) XSpan() M {
 	return e[2] - e[0]
 }
 
 // YSpan is the distance of the Extent in Y
-func (e Extent) YSpan() int64 {
+func (e Extent) YSpan() M {
 	return e[3] - e[1]
 }
 
 // Extent returns back the min and max of the Extent
-func (e Extent) Extent() [4]int64 {
-	return [4]int64{e.MinX(), e.MinY(), e.MaxX(), e.MaxY()}
+func (e Extent) Extent() [4]M {
+	return [4]M{e.MinX(), e.MinY(), e.MaxX(), e.MaxY()}
 }

intgeom/intgeom.go

@@ -1,34 +1,52 @@
 // Package intgeom resembles github.com/go-spatial/geom but uses int64s internally
 // to avoid floating point errors when performing arithmetic with the coords.
+//
 // The idea is that an int64's range (math.MaxInt64) is enough for (most) (earthly) geom operations.
 // See https://www.explainxkcd.com/wiki/index.php/2170:_Coordinate_Precision.
-// Using the last 9 digits as decimals should be enough for identifying
+// Using the last 10 digits as decimals should be enough for identifying
 // the location of a grain of sand (in degrees).
-// That leaves 10 digits for the whole units of measurement in your SRS.
+// More importantly, 10 digits are necessary to minimize the rounding error
+// when the span of a matrix is divided by its size (in pixels) on deeper levels.
+//
+// That leaves 9 digits for the whole units of measurement in your SRS.
 // If that unit is degrees, it's more than enough (a circle only has 360).
 // If that unit is feet, earth's circumference (131 482 560) also still fits.
 // This is not intended to cover everything. go-spatial/geom has much more functionality.
 // You are not required to use this.
 package intgeom
 
-import "github.com/go-spatial/geom/planar"
+import (
+	"fmt"
+	"math"
+	"strconv"
+	"strings"
+
+	"github.com/go-spatial/geom/planar"
+)
 
 const (
 	debug     = false
-	Precision = 1e09
+	Precision = 10
+	Half      = 5000000000
+	One       = 10000000000
 )
 
+// M is short for measure.
+// Used to indicate that a distance or ordinate is saved as an int64 and needs division by Precision (eventually).
+// Shortened for readability in long lines (in other packages it will also be prefixed with intgeom.)
+type M = int64
+
 // ToGeomOrd turns an ordinate represented as an integer back into a floating point
-func ToGeomOrd(o int64) float64 {
+func ToGeomOrd(o M) float64 {
 	if o == 0 {
 		return 0.0
 	}
-	return float64(o) / Precision
+	return float64(o) / math.Pow(10, Precision)
 }
 
 // FromGeomOrd turns a floating point ordinate into a representation by an integer
-func FromGeomOrd(o float64) int64 {
-	return int64(o * Precision)
+func FromGeomOrd(o float64) M {
+	return int64(o * math.Pow(10, Precision))
 }
 
 // SegmentIntersect will find the intersection point (x,y) between two lines if
@@ -41,3 +59,16 @@ func SegmentIntersect(l1, l2 Line) ([2]int64, bool) {
 	intIntersection := [2]int64{FromGeomOrd(intersection[0]), FromGeomOrd(intersection[0])}
 	return intIntersection, intersects
 }
+
+func PrintWithDecimals(o M, n uint) string {
+	s := fmt.Sprintf("%0"+strconv.Itoa(Precision+1)+"d", o)
+	l := len(s)
+	m := s[l-Precision : l]
+	if n < Precision {
+		m = m[0:n]
+	} else {
+		m += strings.Repeat("0", int(n-Precision))
+	}
+	c := s[0 : l-Precision]
+	return c + "." + m
+}

intgeom/line.go

@@ -7,7 +7,7 @@ import (
 )
 
 // Line has exactly two points
-type Line [2][2]int64
+type Line [2][2]M
 
 func (l Line) ToGeomLine() geom.Line {
 	return geom.Line{
@@ -35,10 +35,10 @@ func (l Line) Point1() *Point { return (*Point)(&l[0]) }
 // Point2 returns a new copy of the second point in the line.
 func (l Line) Point2() *Point { return (*Point)(&l[1]) }
 
-func (l Line) Vertices() [][2]int64 { return l[:] }
+func (l Line) Vertices() [][2]M { return l[:] }
 
 // ContainsPoint checks to see if the given pont lines on the line segment. (Including the end points.)
-func (l Line) ContainsPoint(pt [2]int64) bool {
+func (l Line) ContainsPoint(pt [2]M) bool {
 	minx, maxx := l[0][0], l[1][0]
 	if minx > maxx {
 		minx, maxx = maxx, minx

intgeom/point.go

@@ -5,7 +5,7 @@ import (
 )
 
 // Point describes a simple 2D point
-type Point [2]int64
+type Point [2]M
 
 func (p Point) ToGeomPoint() geom.Point {
 	return geom.Point{
@@ -22,18 +22,18 @@ func FromGeomPoint(p geom.Point) Point {
 }
 
 // XY returns an array of 2D coordinates
-func (p Point) XY() [2]int64 {
+func (p Point) XY() [2]M {
 	return p
 }
 
 // SetXY sets a pair of coordinates
-func (p Point) SetXY(xy [2]int64) {
+func (p Point) SetXY(xy [2]M) {
 	p[0] = xy[0]
 	p[1] = xy[1]
 }
 
 // X is the x coordinate of a point in the projection
-func (p Point) X() int64 { return p[0] }
+func (p Point) X() M { return p[0] }
 
 // Y is the y coordinate of a point in the projection
-func (p Point) Y() int64 { return p[1] }
+func (p Point) Y() M { return p[1] }

main.go

@@ -7,8 +7,11 @@ import (
 	"log"
 	"os"
 	"path"
+	"slices"
 	"syscall"
 
+	"github.com/pdok/texel/pointindex"
+
 	"github.com/go-spatial/geom"
 
 	"github.com/carlmjohnson/versioninfo"
@@ -101,6 +104,9 @@ func main() {
 		if err != nil {
 			return err
 		}
+		if err = validateTileMatrixSet(tileMatrixSet, tileMatrixIDs); err != nil {
+			return err
+		}
 
 		_, err = os.Stat(c.String(SOURCE))
 		if os.IsNotExist(err) {
@@ -158,6 +164,19 @@ func main() {
 	}
 }
 
+func validateTileMatrixSet(tms tms20.TileMatrixSet, tileMatrixIDs []tms20.TMID) error {
+	deepestTMID := slices.Max(tileMatrixIDs)
+	stats, deviationInUnits, deviationInPixels, err := pointindex.DeviationStats(tms, deepestTMID)
+	if err != nil {
+		return err
+	}
+	if deviationInPixels >= 1 {
+		log.Printf("warning, (largest) deviation is larger than 1 tile pixel (%f units) on the deepest matrix (%d)\n", deviationInUnits, deepestTMID)
+		log.Println(stats)
+	}
+	return pointindex.IsQuadTree(tms)
+}
+
 func initGPKGTarget(targetPathFmt string, tmID int, overwrite bool, pagesize int) *gpkg.TargetGeopackage {
 	targetPath := fmt.Sprintf(targetPathFmt, tmID)
 	if overwrite {

mathhelp/mathhelp.go

@@ -1,6 +1,13 @@
 package mathhelp
 
-func BetweenInc(f, p, q int64) bool {
+func IBetweenInc(f, p, q int64) bool {
+	if p <= q {
+		return p <= f && f <= q
+	}
+	return q <= f && f <= p
+}
+
+func FBetweenInc(f, p, q float64) bool {
 	if p <= q {
 		return p <= f && f <= q
 	}