diff --git a/intgeom/extent.go b/intgeom/extent.go index 9cda519..f5ee52d 100644 --- a/intgeom/extent.go +++ b/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()} } diff --git a/intgeom/intgeom.go b/intgeom/intgeom.go index 2d733bd..05c4006 100644 --- a/intgeom/intgeom.go +++ b/intgeom/intgeom.go @@ -1,10 +1,14 @@ // 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. @@ -22,13 +26,18 @@ import ( const ( debug = false - Precision = 9 + Precision = 10 + Half = 5000000000 + One = 10000000000 ) -type Ord = int64 +// 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 Ord) float64 { +func ToGeomOrd(o M) float64 { if o == 0 { return 0.0 } @@ -36,7 +45,7 @@ func ToGeomOrd(o Ord) float64 { } // FromGeomOrd turns a floating point ordinate into a representation by an integer -func FromGeomOrd(o float64) Ord { +func FromGeomOrd(o float64) M { return int64(o * math.Pow(10, Precision)) } @@ -51,14 +60,14 @@ func SegmentIntersect(l1, l2 Line) ([2]int64, bool) { return intIntersection, intersects } -func PrintWithDecimals(o Ord, n uint) string { +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 = m + strings.Repeat("0", int(n-Precision)) + m += strings.Repeat("0", int(n-Precision)) } c := s[0 : l-Precision] return c + "." + m diff --git a/intgeom/line.go b/intgeom/line.go index 46857e9..395a91a 100644 --- a/intgeom/line.go +++ b/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 diff --git a/intgeom/point.go b/intgeom/point.go index f14f934..20ab826 100644 --- a/intgeom/point.go +++ b/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] } diff --git a/main.go b/main.go index 692f638..6fd4c95 100644 --- a/main.go +++ b/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 nil +} + func initGPKGTarget(targetPathFmt string, tmID int, overwrite bool, pagesize int) *gpkg.TargetGeopackage { targetPath := fmt.Sprintf(targetPathFmt, tmID) if overwrite { diff --git a/pointindex/pointindex.go b/pointindex/pointindex.go index 4bf71bd..4b54cf2 100644 --- a/pointindex/pointindex.go +++ b/pointindex/pointindex.go @@ -1,10 +1,13 @@ package pointindex import ( + "errors" "fmt" + "golang.org/x/exp/maps" "io" "math" "slices" + "strconv" "github.com/pdok/texel/mathhelp" "github.com/pdok/texel/morton" @@ -30,8 +33,6 @@ const ( // bottomright = bottom | right // 0b01 // topleft = top | left // 0b10 // topright = top | right // 0b11 - intHalf = 500000000 - intOne = 1000000000 ) type Quadrant struct { @@ -62,9 +63,10 @@ type Quadrant struct { // inc type PointIndex struct { Quadrant - maxDepth Level - // Number of quadrants (in one direction) on the deepest level (= 2 ^ maxDepth) + deepestLevel Level + // Number of quadrants (in one direction) on the deepest level (= 2 ^ deepestLevel) deepestSize uint + deepestRes intgeom.M quadrants map[Level]map[morton.Z]Quadrant hitOnce map[Level]map[intgeom.Point][]int hitMultiple map[Level]map[intgeom.Point][]int @@ -74,10 +76,12 @@ type Level = uint type Q = int // quadrant index (0, 1, 2 or 3) func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID) *PointIndex { - // TODO ensure that the tile matrix set is actually a quad tree, starting at 0. just assuming now. + if err := isQuadTree(tileMatrixSet); err != nil { + panic(err) + } rootTM := tileMatrixSet.TileMatrices[0] levelDiff := uint(math.Log2(float64(rootTM.TileWidth))) + uint(math.Log2(float64(VectorTileInternalPixelResolution))) - maxDepth := uint(deepestTMID) + levelDiff + deepestLevel := uint(deepestTMID) + levelDiff bottomLeft, topRight, err := tileMatrixSet.MatrixBoundingBox(0) if err != nil { panic(fmt.Errorf(`could not make PointIndex from TileMatrixSet %v: %w'`, tileMatrixSet.ID, err)) @@ -85,18 +89,20 @@ func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID intBottomLeft := intgeom.FromGeomPoint(bottomLeft) intTopRight := intgeom.FromGeomPoint(topRight) intExtent := intgeom.Extent{intBottomLeft.X(), intBottomLeft.Y(), intTopRight.X(), intTopRight.Y()} + deepestSize := mathhelp.Pow2(deepestLevel) ix := PointIndex{ Quadrant: Quadrant{ intExtent: intExtent, z: 0, }, - maxDepth: maxDepth, - deepestSize: mathhelp.Pow2(maxDepth), - quadrants: make(map[Level]map[morton.Z]Quadrant, maxDepth+1), - hitOnce: make(map[uint]map[intgeom.Point][]int), - hitMultiple: make(map[uint]map[intgeom.Point][]int), + deepestLevel: deepestLevel, + deepestSize: deepestSize, + deepestRes: intExtent.XSpan() / int64(deepestSize), + quadrants: make(map[Level]map[morton.Z]Quadrant, deepestLevel+1), + hitOnce: make(map[uint]map[intgeom.Point][]int), + hitMultiple: make(map[uint]map[intgeom.Point][]int), } - _, ix.intCentroid = getQuadrantExtentAndCentroid(0, 0, 0, intExtent) + _, ix.intCentroid = ix.getQuadrantExtentAndCentroid(0, 0, 0, intExtent) return &ix } @@ -109,7 +115,7 @@ func (ix *PointIndex) InsertPolygon(polygon geom.Polygon) { pointsCount += len(ring) } var level uint - for level = 0; level <= ix.maxDepth; level++ { + for level = 0; level <= ix.deepestLevel; level++ { if ix.quadrants[level] == nil { ix.quadrants[level] = make(map[morton.Z]Quadrant, pointsCount) // TODO smaller for the shallower levels } @@ -125,10 +131,9 @@ func (ix *PointIndex) InsertPolygon(polygon geom.Polygon) { // InsertPoint inserts a Point by its absolute coord func (ix *PointIndex) InsertPoint(point geom.Point) { - intDeepestRes := ix.intExtent.XSpan() / int64(ix.deepestSize) intPoint := intgeom.FromGeomPoint(point) - deepestX := int((intPoint.X() - ix.intExtent.MinX()) / intDeepestRes) - deepestY := int((intPoint.Y() - ix.intExtent.MinY()) / intDeepestRes) + deepestX := int((intPoint.X() - ix.intExtent.MinX()) / ix.deepestRes) + deepestY := int((intPoint.Y() - ix.intExtent.MinY()) / ix.deepestRes) ix.InsertCoord(deepestX, deepestY) } @@ -144,14 +149,14 @@ func (ix *PointIndex) InsertCoord(deepestX int, deepestY int) { // insertCoord adds a point into this pc, assuming the point is inside its extent func (ix *PointIndex) insertCoord(deepestX int, deepestY int) { var l Level - for l = 0; l <= ix.maxDepth; l++ { - x := uint(deepestX) / mathhelp.Pow2(ix.maxDepth-l) - y := uint(deepestY) / mathhelp.Pow2(ix.maxDepth-l) + for l = 0; l <= ix.deepestLevel; l++ { + x := uint(deepestX) / mathhelp.Pow2(ix.deepestLevel-l) + y := uint(deepestY) / mathhelp.Pow2(ix.deepestLevel-l) z := morton.MustToZ(x, y) if ix.quadrants[l] == nil { // probably already initialized by InsertPolygon ix.quadrants[l] = make(map[morton.Z]Quadrant) } - extent, centroid := getQuadrantExtentAndCentroid(l, x, y, ix.intExtent) + extent, centroid := ix.getQuadrantExtentAndCentroid(l, x, y, ix.intExtent) ix.quadrants[l][z] = Quadrant{ z: z, intExtent: extent, @@ -160,8 +165,8 @@ func (ix *PointIndex) insertCoord(deepestX int, deepestY int) { } } -func getQuadrantExtentAndCentroid(level Level, x, y uint, intRootExtent intgeom.Extent) (intgeom.Extent, intgeom.Point) { - intQuadrantSpan := intRootExtent.XSpan() / int64(mathhelp.Pow2(level)) +func (ix *PointIndex) getQuadrantExtentAndCentroid(level Level, x, y uint, intRootExtent intgeom.Extent) (intgeom.Extent, intgeom.Point) { + intQuadrantSpan := int64(mathhelp.Pow2(ix.deepestLevel-level)) * ix.deepestRes intMinX := intRootExtent.MinX() intMinY := intRootExtent.MinY() intExtent := intgeom.Extent{ @@ -218,7 +223,7 @@ func (ix *PointIndex) snapClosestPoints(intLine intgeom.Line, levelMap map[Level } var level Level - for level = 1; level <= ix.maxDepth; level++ { + for level = 1; level <= ix.deepestLevel; level++ { quadrantsIntersected := make([]Quadrant, 0, 10) // TODO good estimate of expected count, based on line length / quadrant span * geom points? for _, parent := range parents { quadrantZs := getQuadrantZs(parent.z) @@ -483,10 +488,88 @@ func (ix *PointIndex) ToWkt(writer io.Writer) { for _, quadrant := range quadrants { _ = wkt.Encode(writer, quadrant.intExtent.ToGeomExtent()) _, _ = fmt.Fprintf(writer, "\n") - if level == ix.maxDepth { + if level == ix.deepestLevel { _ = wkt.Encode(writer, quadrant.intCentroid.ToGeomPoint()) _, _ = fmt.Fprintf(writer, "\n") } } } } + +func isQuadTree(tms tms20.TileMatrixSet) error { + var previousTMID int + var previousTM *tms20.TileMatrix + tmIDs := maps.Keys(tms.TileMatrices) + slices.Sort(tmIDs) + for _, tmID := range tmIDs { + tm := tms.TileMatrices[tmID] + if tm.MatrixHeight != tm.MatrixWidth { + return errors.New("tile matrix height should be same as width: " + tm.ID) + } + if tm.TileHeight != tm.TileWidth { + return errors.New("tiles should be square: " + tm.ID) + } + tmIDStringToInt, err := strconv.Atoi(tm.ID) + if err != nil { + return err + } + if tmIDStringToInt != tmID { + return errors.New("tile matrix ID should string representation of its index in the array: " + tm.ID) + } + if len(tm.VariableMatrixWidths) != 0 { + return errors.New("variable matrix widths are not supported: " + tm.ID) + } + if previousTM != nil { + if tmID != previousTMID+1 { + return errors.New("tile matrix IDs should be a range with step 1 starting with 0") + } + if *tm.PointOfOrigin != *previousTM.PointOfOrigin { + return errors.New("tile matrixes should have the same point of origin: " + tm.ID) + } + if tm.CornerOfOrigin != previousTM.CornerOfOrigin { + return errors.New("tile matrixes should have the same corner of origin: " + tm.ID) + } + } + + previousTMID = tmID + previousTM = &tm + } + return nil +} + +// DeviationStats calcs some stats to show the result of using ints internally +// if the res (span/size) is not "round" (with intgeom.Precision), a deviation will arise +// +//nolint:nakedret +func DeviationStats(tms tms20.TileMatrixSet, deepestTMID tms20.TMID) (stats string, deviationInUnits, deviationInPixels float64, err error) { + bottomLeft, topRight, err := tms.MatrixBoundingBox(0) + if err != nil { + return + } + ix := FromTileMatrixSet(tms, deepestTMID) + p := uint(intgeom.Precision + 1) + ps := strconv.Itoa(int(p)) + stats += fmt.Sprintf("deepest level: %d\n", ix.deepestLevel) + stats += fmt.Sprintf("deepest pixels count on one axis: %d\n", ix.deepestSize) + stats += fmt.Sprintf("float minx, miny: %."+ps+"f, %."+ps+"f\n", bottomLeft.X(), bottomLeft.Y()) + stats += fmt.Sprintf("float maxx, maxy: %."+ps+"f, %."+ps+"f\n", topRight.X(), topRight.Y()) + floatSpanX := topRight.X() - bottomLeft.X() + stats += fmt.Sprintf("float span X: %."+ps+"f\n", floatSpanX) + stats += fmt.Sprintf("int64 span X: %s\n", intgeom.PrintWithDecimals(ix.intExtent.XSpan(), p)) + floatRes := floatSpanX / float64(ix.deepestSize) + stats += fmt.Sprintf("float reso: %."+ps+"f\n", floatRes) + intRes := ix.intExtent.XSpan() / int64(ix.deepestSize) + stats += fmt.Sprintf("int64 reso: %s\n", intgeom.PrintWithDecimals(intRes, p)) + + floatRecalcMaxX := floatRes * float64(ix.deepestSize) + stats += fmt.Sprintf("float recalc maxX: %."+ps+"f\n", floatRecalcMaxX) + intRecalcMaxX := intgeom.ToGeomOrd(intRes * int64(ix.deepestSize)) + stats += fmt.Sprintf("int64 recalc maxX: %."+ps+"f\n", intRecalcMaxX) + + deviationInUnits = floatRecalcMaxX - intRecalcMaxX + deviationInPixels = deviationInUnits / floatRes + stats += fmt.Sprintf("deviation (in units) maxX : %."+ps+"f\n", deviationInUnits) + stats += fmt.Sprintf("deviation (in pixels) maxX: %."+ps+"f\n", deviationInPixels) + + return +} diff --git a/pointindex/pointindex_test.go b/pointindex/pointindex_test.go index 81bdc12..fc3eb44 100644 --- a/pointindex/pointindex_test.go +++ b/pointindex/pointindex_test.go @@ -71,7 +71,7 @@ func TestPointIndex_containsPoint(t *testing.T) { want: false, }, } - extent := intgeom.Extent{0, 0, intOne, intOne} + extent := intgeom.Extent{0, 0, intgeom.One, intgeom.One} for _, tt := range tests { t.Run(tt.name, func(t *testing.T) { got := containsPoint(intgeom.FromGeomPoint(tt.pt), extent) @@ -94,10 +94,10 @@ func TestPointIndex_getQuadrantExtentAndCentroid(t *testing.T) { }{ { name: "simple", - intRootExtent: intgeom.Extent{0, 0, intOne, intOne}, + intRootExtent: intgeom.Extent{0, 0, intgeom.One, intgeom.One}, want: wants{ - extent: intgeom.Extent{0, 0, intOne, intOne}, - centroid: intgeom.Point{intHalf, intHalf}, + extent: intgeom.Extent{0, 0, intgeom.One, intgeom.One}, + centroid: intgeom.Point{intgeom.Half, intgeom.Half}, }, }, { @@ -111,7 +111,12 @@ func TestPointIndex_getQuadrantExtentAndCentroid(t *testing.T) { } for _, tt := range tests { t.Run(tt.name, func(t *testing.T) { - extent, centroid := getQuadrantExtentAndCentroid(0, 0, 0, tt.intRootExtent) + ix := PointIndex{ + deepestLevel: 0, + deepestSize: 1, + deepestRes: tt.intRootExtent.XSpan() / 1, + } + extent, centroid := ix.getQuadrantExtentAndCentroid(0, 0, 0, tt.intRootExtent) if !assert.EqualValues(t, tt.want.extent, extent) { t.Errorf("getQuadrantExtentAndCentroid() = %v, want %v", extent, tt.want.extent) } @@ -138,8 +143,9 @@ func TestPointIndex_InsertPoint(t *testing.T) { intExtent: intgeom.FromGeomExtent(geom.Extent{0.0, 0.0, 1.0, 1.0}), intCentroid: intgeom.FromGeomPoint(geom.Point{0.5, 0.5}), }, - maxDepth: 0, - deepestSize: mathhelp.Pow2(0), + deepestLevel: 0, + deepestSize: mathhelp.Pow2(0), + deepestRes: intgeom.FromGeomOrd(1.0) / intgeom.M(mathhelp.Pow2(0)), quadrants: map[Level]map[morton.Z]Quadrant{0: {0: Quadrant{ intExtent: intgeom.FromGeomExtent(geom.Extent{0.0, 0.0, 1.0, 1.0}), intCentroid: intgeom.FromGeomPoint(geom.Point{0.5, 0.5}), @@ -155,8 +161,9 @@ func TestPointIndex_InsertPoint(t *testing.T) { intExtent: intgeom.FromGeomExtent(geom.Extent{0.0, 0.0, 1.0, 1.0}), intCentroid: intgeom.FromGeomPoint(geom.Point{0.5, 0.5}), }, - maxDepth: 1, - deepestSize: mathhelp.Pow2(1), + deepestLevel: 1, + deepestSize: mathhelp.Pow2(1), + deepestRes: intgeom.FromGeomOrd(1.0) / intgeom.M(mathhelp.Pow2(1)), quadrants: map[Level]map[morton.Z]Quadrant{ 0: {0: Quadrant{ intExtent: intgeom.FromGeomExtent(geom.Extent{0.0, 0.0, 1.0, 1.0}), @@ -179,8 +186,9 @@ func TestPointIndex_InsertPoint(t *testing.T) { intExtent: intgeom.FromGeomExtent(geom.Extent{0.0, 0.0, 4.0, 4.0}), intCentroid: intgeom.FromGeomPoint(geom.Point{2.0, 2.0}), }, - maxDepth: 3, - deepestSize: mathhelp.Pow2(3), + deepestLevel: 3, + deepestSize: mathhelp.Pow2(3), + deepestRes: intgeom.FromGeomOrd(4.0) / intgeom.M(mathhelp.Pow2(3)), quadrants: map[Level]map[morton.Z]Quadrant{ 0: {0: Quadrant{ z: 0, @@ -214,8 +222,9 @@ func TestPointIndex_InsertPoint(t *testing.T) { intExtent: intgeom.FromGeomExtent(geom.Extent{0.0, 0.0, 16.0, 16.0}), intCentroid: intgeom.FromGeomPoint(geom.Point{8.0, 8.0}), }, - maxDepth: 5, - deepestSize: mathhelp.Pow2(5), + deepestLevel: 5, + deepestSize: mathhelp.Pow2(5), + deepestRes: intgeom.FromGeomOrd(16.0) / intgeom.M(mathhelp.Pow2(5)), quadrants: map[Level]map[morton.Z]Quadrant{ 0: {0: Quadrant{ z: 0, @@ -314,8 +323,8 @@ func TestPointIndex_InsertPoint_Deepest(t *testing.T) { ix := FromTileMatrixSet(tms, tt.tmID) ix.InsertPoint(tt.point) - assert.Equal(t, 1, len(ix.quadrants[ix.maxDepth])) - for z, quadrant := range ix.quadrants[ix.maxDepth] { + assert.Equal(t, 1, len(ix.quadrants[ix.deepestLevel])) + for z, quadrant := range ix.quadrants[ix.deepestLevel] { assert.EqualValues(t, tt.want, quadrant) assert.EqualValues(t, tt.want.z, z) } @@ -455,7 +464,7 @@ func TestPointIndex_SnapClosestPoints(t *testing.T) { ix.InsertPolygon(poly) levels := tt.levels if levels == nil { - levels = []Level{ix.maxDepth} + levels = []Level{ix.deepestLevel} } got := ix.SnapClosestPoints(tt.line, mapslicehelp.AsKeys(levels), tt.ringID) if !assert.EqualValues(t, tt.want, got) { @@ -499,19 +508,22 @@ func TestPointIndex_lineIntersects(t *testing.T) { } } -func newSimplePointIndex(maxDepth Level, cellSize float64) *PointIndex { - span := cellSize * float64(mathhelp.Pow2(maxDepth)) +func newSimplePointIndex(deepestLevel Level, cellSize float64) *PointIndex { + deepestSize := mathhelp.Pow2(deepestLevel) + span := cellSize * float64(deepestSize) + intExtent := intgeom.Extent{0.0, 0.0, intgeom.FromGeomOrd(span), intgeom.FromGeomOrd(span)} ix := PointIndex{ Quadrant: Quadrant{ - intExtent: intgeom.Extent{0.0, 0.0, intgeom.FromGeomOrd(span), intgeom.FromGeomOrd(span)}, - }, - maxDepth: maxDepth, - deepestSize: mathhelp.Pow2(maxDepth), - quadrants: make(map[Level]map[morton.Z]Quadrant, maxDepth+1), - hitOnce: make(map[morton.Z]map[intgeom.Point][]int, 0), - hitMultiple: make(map[morton.Z]map[intgeom.Point][]int, 0), + intExtent: intExtent, + }, + deepestLevel: deepestLevel, + deepestSize: deepestSize, + deepestRes: intExtent.XSpan() / int64(deepestSize), + quadrants: make(map[Level]map[morton.Z]Quadrant, deepestLevel+1), + hitOnce: make(map[morton.Z]map[intgeom.Point][]int, 0), + hitMultiple: make(map[morton.Z]map[intgeom.Point][]int, 0), } - _, ix.intCentroid = getQuadrantExtentAndCentroid(0, 0, 0, ix.intExtent) + _, ix.intCentroid = ix.getQuadrantExtentAndCentroid(0, 0, 0, ix.intExtent) return &ix } diff --git a/snap/snap_test.go b/snap/snap_test.go index 88f0692..2aadce9 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -915,18 +915,18 @@ func Test_dedupeInnersOuters(t *testing.T) { } // newSimpleTileMatrixSet creates a tms for snap testing purposes -// the effective quadrant amount (one axis) on the deepest level will be 2^maxDepth * 16 (vt internal pixel res) +// the effective quadrant amount (one axis) on the deepest level will be 2^deepestLevel * 16 (vt internal pixel res) // the effective quadrant size (one axis) on the deepest level will be cellSize / 16 (vt internal pixel res) -func newSimpleTileMatrixSet(maxDepth pointindex.Level, cellSize float64) tms20.TileMatrixSet { +func newSimpleTileMatrixSet(deepestTMID pointindex.Level, cellSize float64) tms20.TileMatrixSet { zeroZero := tms20.TwoDPoint([2]float64{0.0, 0.0}) tms := tms20.TileMatrixSet{ CRS: fakeCRS{}, OrderedAxes: []string{"X", "Y"}, - TileMatrices: make(map[tms20.TMID]tms20.TileMatrix, maxDepth+1), + TileMatrices: make(map[tms20.TMID]tms20.TileMatrix, deepestTMID+1), } - for tmID := 0; tmID <= int(maxDepth); tmID++ { + for tmID := 0; tmID <= int(deepestTMID); tmID++ { // (only values from the root tm are used, for the rest it is assumed to follow quad matrix rules) - tmCellSize := cellSize * float64(mathhelp.Pow2(maxDepth-uint(tmID))) + tmCellSize := cellSize * float64(mathhelp.Pow2(deepestTMID-uint(tmID))) tms.TileMatrices[tmID] = tms20.TileMatrix{ ID: "0", ScaleDenominator: tmCellSize / tms20.StandardizedRenderingPixelSize, @@ -942,7 +942,6 @@ func newSimpleTileMatrixSet(maxDepth pointindex.Level, cellSize float64) tms20.T return tms } -//nolint:unparam func loadEmbeddedTileMatrixSet(t *testing.T, tmsID string) tms20.TileMatrixSet { tms, err := tms20.LoadEmbeddedTileMatrixSet(tmsID) require.NoError(t, err)