go-spatial/geom

Geometry interfaces to help drive interoperability within the Go geospatial community

Packages

package geom

Geometry interfaces to help drive interoperability within the Go geospatial community. This package focuses on 2D geometries.

package encoding

Encoding package describes a few interfaces and common errors for various sub packages implementing different encoders.

Dependencies

Dependencies are managed through go mod with the execption of one package:

• https://github.com/dhconnelly/rtreego BSD 3 License We are keeping this internal, so that we can build an rtree implementation that uses geom types as it's base, but is build ontop of this.

geom's Issues

Polygon rendering issue. Not sure how to describe [screenshots]

Tile ZXY: /4/14/6

The first screenshot shows a broken tile. At first glance, it looks like a winding order issue but if you look at the second screenshot the area that is a cutout shows a feature when hovered. When the inspector encounters a cutout no feature information should be displayed.

The following data is the entire planet which is used at zoom 4. This is Natural Earth data, which has already been simplified (no simplification happening in tegola)

data-1573683452170.csv.zip

Add a testgeometry directory

When using this package in other directories it would be nice to have some pre-made geometries to use for testing compiler and fatal runtime errors. I propose to put them in a subpackage of the geom repo, something like gtest.

Add members to the geometry interface

The geometry interface is not very helpful since any type can match to it and we should use the compiler enforcement when possible.

possibilities:

type Geometry interface {
  Points() []geom.Point
  Clone() geom.Geometry
  WalkPoints(func(*geom.Point) (err error, ok bool)) error
  String() string
}

Point within Polygon?

I dont see a function for finding whether a point is within a polygon. Am I missing it, or can I request it be added to the API.

Thanks

wkt_decode: readWhitespace does not skip all white space.

readWhitespace assumes all text is ASCII. The function should decode runes correctly and then check to see if they are spaces.

Can not test runes this way:

https://github.com/go-spatial/geom/blob/master/encoding/wkt/wkt_decode.go#L55-L57

You need to read all the bytes in the run before testing it.

Here is a test case with the failures:

func TestReadWhitespace(t *testing.T) {
	type tcase struct {
		in    string
		match bool
	}
	fn := func(tc tcase) func(*testing.T) {
		decoder := NewDecoder(strings.NewReader(tc.in + "."))

		return func(t *testing.T) {
			match, err := decoder.readWhitespace()
			if err != nil {
				t.Errorf("error, expected nil got: %v", err)
				return
			}
			if match != tc.match {
				t.Errorf("match, expected %t, got %t", tc.match, match)
			}
		}
	}

	tests := map[string]tcase{
		"regular0":     tcase{in: " ", match: true},
		"tab":          tcase{in: "\t", match: true},
		"newline":      tcase{in: "\n", match: true},
		"return":       tcase{in: "\r", match: true},
		"regular":      tcase{in: "\u0020", match: true},
		"nobreak":      tcase{in: "\u00A0", match: true},
		"ogham mark":   tcase{in: "\u1680", match: true},
		"en quad":      tcase{in: "\u2000", match: true},
		"em quad":      tcase{in: "\u2001", match: true},
		"en":           tcase{in: "\u2002", match: true},
		"em":           tcase{in: "\u2003", match: true},
		"three-per-em": tcase{in: "\u2004", match: true},
		"four-per-em":  tcase{in: "\u2005", match: true},
		"six-per-em":   tcase{in: "\u2006", match: true},
		"figure":       tcase{in: "\u2007", match: true},
		"punctuation":  tcase{in: "\u2007", match: true},
		"thin":         tcase{in: "\u2009", match: true},
		"hair":         tcase{in: "\u200A", match: true},
		"nnbsp":        tcase{in: "\u202F", match: true},
		"mmsp":         tcase{in: "\u205F", match: true},
		"ideographic":  tcase{in: "\u3000", match: true},
	}
	for key, test := range tests {
		t.Run(key, fn(test))
	}

}

WKT Decoder

WKT Encoder is complete. Now we need a WKT Decoder

Naming of mvt.ScaleGeo

I also came across some code in mvt/scale.go that does not do strictly scaling (it also removes points, which is slightly misleading), and @ARolek and I came up with the idea of renaming it PrepareGeo or something similar. This could do some sort of scaling, de-duplication, and standardize how we handle something like the last point of polygon.

Another thought is that scaling could be part of the validation routine as a whole.

Note that this issue might be moved to the go-spatial/geom repo.

[geom] dropping support for interface types

Support for things like Pointer, LineStringer, etc. Should be dropped from the general code base. They were added for interoperability, but support has been sparse and adds an extra cognitive load in switch x.(type) using code.

For interoperability, I propose a convert or interop package with contains the interfaces and and has a ToGeom function. A readme would also have an example on how to write a FromGeom function, and perhaps a FromGeomer interface which contains the necessary functions to convert from the geom types.

Implement constraints for the triangulation

Currently, constraints are not working.

wkt: behavior for pointless geometries and empty

Currently, the encoder will clean up "pointless" (hehe) geometries and count them as empty.
For example, from the wkt tests:

{
    Geom: &geom.MultiPolygon{{{}}, {{}}},
    Rep:  "MULTIPOLYGON EMPTY",
},

At the moment, I'm cleaning up the wkt parser and my, simpler, implementation does not do this (yet?). I've been using the postgis the implementations of a wkt {en,de}coder and it seems that it is inconsistent with itself. I believe the fidelity of the structure of the geometries should be kept, but it seems the convention is to remove pointless geometries.

postgis output:

$ psql -d natural_earth -t -c "select ST_AsText(ST_AsBinary('geometrycollection ( linestring empty )'::geometry)::geometry)"
 GEOMETRYCOLLECTION EMPTY
$ psql -d natural_earth -t -c "select ST_AsText(ST_GeomFromText('geometrycollection (linestring empty)'))"   
 GEOMETRYCOLLECTION(LINESTRING EMPTY)

Check to see if the constraints provided for triangulation are colinear.

If the constraints provided for triangulation are colinear we should combine them into a single constraint.

mvt: Remove tegola.Tile from prepare function

While writing the mvt docs, I realized I did not remove tegola.Tile from this function

geom/encoding/mvt/prepare.go

Line 11 in b1d0cb8

func PrepareGeo(geo tegola.Geometry, tile *tegola.Tile) geom.Geometry {

[makevalid] Cutout ring has extra point, causing extra line to be drawn:

This is using the natural_earth water dataset.
Location of error:
http://localhost:8080/?debug=true#13.13/52.75645/173.95104

Polygon In Question
ClipArea: slippy.NewTile(13,8054,2677).Extent3857().ExpandBy(64.0)

The issue is that one of the "bridges" is remaining in the cutout sliver.

This is an issue in the PolygonForRing function

[triangulation] Failed to insert point

Another failed to insert point from scoutred data set
main_points.wkt.txt

wkt encodes float with lowercase 'e' scientific format, but decode() expects uppercase 'E'

In wkt_encode.go, default format is 'g'. 0.0000000000001 will be encoded as '1e-12'

// defaults of strict = false, precision = 10, and fmt = 'g'.
func NewDefaultEncoder(w io.Writer) Encoder {
	return NewEncoder(w, false, 10, 'g')
}

wkt_decode.go expects 'E':

	isNumeric := func(b byte) bool {
		return (b >= '0' && b <= '9') ||
			b == '-' ||
			b == '.' ||
			// b == ',' || // technically part of the spec,
			// but even postgis does not support it
			b == 'E'
	}
...

mvt: handling south- and west-positive projection system

mvt.PrepareGeometry (and mvt as whole) was written to be srid agnostic, but the follwing lines would fail in the edge case that the projection had axis oriented differently from the common web mercator and lat/longs, such as:

• epsg:2053
• epsg:22385
• an antiquated coordinate for Mars is west-positive

These types of coordinate systems are also refered to as Lo coordinate systems.

The following lines are relevant the mvt code:

geom/encoding/mvt/prepare.go

Lines 66 to 67 in f76f5a1

	px := int64((g.X() - tile.MinX()) / tile.XSpan() * pixelExtent)
	py := int64((tile.MaxY() - g.Y()) / tile.YSpan() * pixelExtent)

However, I think the change would need a new or re-designed types. One thing that comes to mind is tweaking the Extent elements to be left, top, right, bottom instead of min/max X/Y. This would not require any changes under our current assumptions in tegola of only using epsg 3857 and 4326.

Here are some relevant links to discussions on this topic:

https://community.esri.com/thread/185667-why-is-the-projection-upside-down
https://isis.astrogeology.usgs.gov/IsisSupport/index.php?topic=370.0
http://www.ngi.gov.za/index.php/technical-information/geodesy-and-gps/datum-s-and-coordinate-systems

[makevalid] splitIntersectingLines is dropping too many segments.

In the gdey_delaunay branch.
In the makevalid package the function splitIntersectingLines is dropping too many lines after it has
found the intersecting lines.

The following image: The purple lines are the lines from the original geometry. The gray lines are the segments of the geometry within the clip box and the segments of the clipbox.
The dashed lines are the lines the current splitIntersectingLines function is generating.
The red lines are the missing lines.

[encoding] Add MVT encoding package

MVT encoding is already in tegola, so this issue is mainly around porting the encoding package from the tegola codebase, adopting the core geom types.

Unit tests are failing for Trianglulator

It looks like the unit tests are failing for Trianglulator around the inserting of an index. I have gone down the stack trace a ways and they all seem to lead to the same part of code.

Right now it looks to be around the method insertConstraints().

It seems to affect test cases 2,3 in makeValid and the last one in Trianglulator.

Will try and work up a fix as I figure out how this code works, if you have any ideas that would be great as well.

Reverse Geocoding

howto build reverse geocoding with many polygon and polylinestring as source of data, the part i want to build is:

1. indexing many polygon, and query like ST_Contains(POINT, GeometryField)
2. indexing multiple PolylineString and query like ST_Intersects(Box,GeometryField)

TIA

Helper function for retrieving the underlying type of geom.Geometry

A helper function that would accept the geom.Geometry container type and then return either a string or a GeomType. This is useful when writing an error message like:

fmt.Printf("expected %t, got %t", geom.LineString{}, goem.Type(geom.Geometry))

Potential function signature

func geom.Type(geom.Geometry) GeomType {
  // switch on the type and return it
}

Duplicate functionallity: geom.IsColiner in circle.go

I think one can use windingorder.PointsFor(...).IsColinear() to test to see if the three points are colinear.

I don't think we need geom.IsColinear anymore.

Add support for extended WKB encoding and PostGIS EWKB.

Here is postgis documentation on their extention to wkb to support ZM and srid values.
https://github.com/postgis/postgis/blob/2.1.0/doc/ZMSgeoms.txt

Wiki page on WKB with other consts.

https://en.wikipedia.org/wiki/Well-known_text#Well-known_binary

do this suport wgs84

I want to get line distance and it is wgs84

Build flags for debugger

Currently the debugger cannot be built on macos, and during runtime, it is dependent on spacialite. Adding a build flag would allow this package's test to be used without spacialite.

slippy.RangeFamilyAt doesn't return self or parents

Problem

RangeFamilyAt was rewritten in PR #99 and its behavior fundamentally changed. Previously, it returned:

1. The tile itself, plus
2. If zoom is less than the tile's zoom, all tiles ancestors at the specified zoom,
3. If zoom is greater than the tile's zoom, all tile's descendants at the specified zoom.

After PR #99 it returns only (3). In addition, floating point error introduced by FromNative and ToNative cause the function's output to be unpredictable and therefore difficult to test (see #127).

PR #99 commented out the failing test: d04b99a#diff-e7c11f2955e26a4cf0f286e7b5eb5da15c1e475bf3fb386ad134449680f31354R56

Solution

Unknown. I don't know the background of this function, but Tegola's tests expect it to function the old way, not the new way.

Inserting point into triangulation failure

Here is another inserting point into triangulation failure.

https://gist.github.com/gdey/6f4bf7a564cf82c63acb3893b5789e1c

Why are the UTM coordinates rounded?

Hello,

I'm using the UTM coordinates package of geom and I found out that the fromLngLat method rounds the output UTM coordinates to the meter. My usage needs a resolution down to the centimeter, so it's not enough.
Is there any particular reason why geom is doing this rounding instead of leaving it to the user?

Thanks!
Ben

Fuzzing wkb encoding causes out of memory errors #384

Ref: go-spatial/tegola#384

Fuzzing geom/encoding/wkb (as per #53) reveals a potential bug with several decode functions. I have attached the relevant test output - 010c1c076fce1acae227d578d230bd713dd44528.txt

https://github.com/go-spatial/tegola/blob/2d90e68580efc38a85e3b6c4ad020a80bf66095a/geom/encoding/wkb/internal/decode/decode.go#L155

In this case, the out of memory error occurs when the Collection function reads the number of items num from the WKB byte stream. In this case it is because of a fuzzed value. When thinking of how to sanitize the input I initially thought it might be sufficient to check it was within a sane range. As far as I can tell, however, there is no documented upper limit for the number of Points in a given geometry (Polygons, Collections etc) in WKB other than math.MaxUint32 (4294967295).

This also affects the other geometry decoding functions that use the value returned from the WKB byte stream to do slice allocation:

https://github.com/go-spatial/tegola/blob/2d90e68580efc38a85e3b6c4ad020a80bf66095a/geom/encoding/wkb/internal/decode/decode.go#L134

https://github.com/go-spatial/tegola/blob/2d90e68580efc38a85e3b6c4ad020a80bf66095a/geom/encoding/wkb/internal/decode/decode.go#L120

https://github.com/go-spatial/tegola/blob/2d90e68580efc38a85e3b6c4ad020a80bf66095a/geom/encoding/wkb/internal/decode/decode.go#L99

https://github.com/go-spatial/tegola/blob/2d90e68580efc38a85e3b6c4ad020a80bf66095a/geom/encoding/wkb/internal/decode/decode.go#L78

https://github.com/go-spatial/tegola/blob/2d90e68580efc38a85e3b6c4ad020a80bf66095a/geom/encoding/wkb/internal/decode/decode.go#L64

https://github.com/go-spatial/tegola/blob/2d90e68580efc38a85e3b6c4ad020a80bf66095a/geom/encoding/wkb/internal/decode/decode.go#L41

geojson.Geometry.UnmarshalJSON() panics if JSON doesn't have the right format

I had an issue with a JSON file where the field "type" was missing.

The function func (geo *Geometry) UnmarshalJSON(b []byte) error panics when it tries to dereference geojsonMap["type"].

I think it will be the same if "geometries" or "coordinates" is missing.

The function must check that key exists before dereferencing geojsonMap[key].

Thanks !

types like `type LineString = []Point`

Thank you starting this project. I agree that there should be interoperability within the Go geospatial community. I'm not really sure what that looks like, but I think a set of defined types that everyone uses is a good start.

I am suggesting/wondering why you didn't choose to do something like:

type Point [2]float64
type MultiPoint []Point

type LineString []Point
type MultiLineString []LineString

type Ring LineString
type Polygon []Ring
type MultiPolygon []Polygon

The reason I like this is because you can do something like this with no type assertion and you still get the safety of knowing it's a Ring vs a MultiPoint vs a LineString:

poly = Polygon{....}
outerRingLength = lonlat.Length(poly[0])

or

mls = MultiLineString{ls}
mls = append(mls, otherLineString)

It feels really natural to me because you can use all the build ins like len, append and [:]. If you want to explicitly convert a LineString to a Ring you can just do it Ring(ls)

I have an attempt in paulmach/orb that went one step further to separate planar/projected types from geo/earth types so that if you did line.Length() it would use the correct math. But that was a pain to work with. One you were forced to copy on project. Two, without generics the code to "crop to bounding box" (which is basically the same in both domains) required copy pasting or slowing it all down by wrapping it in a weird/unnecessary interface.

The demo/version I like now is at paulmach/geo. It has one set of geo.* types and then planar and lonlat sub packages for specific operations that are different, e.g. Area, Length, DistanceFrom, etc. If you have two points p1 and p2 the distances the code would look like:

planar.Distance(p1, p2)
lonlat.Distance(p1, p2)

While it's not type safe, it still reads really well, and building stuff like convexhull doesn't require copied code for the different types. The result is in the same space as the input.

Anyway, this is stuff I've been thinking a lot about lately and from this repo and twitter it seem like you guys have to. Interest to know if you have any thoughts.

slippy.PixelsToNative ignores pixels argument

Problem

PixelsToNative should convert the accepted pixels number into projected/native units. Instead, it doesn't use the variable at all and merely returns the ratio of native to projected.

Solution

Refactor PixelsToNative to either:

1. Convert pixels, or
2. Not accept pixels (in which case, the function should be renamed to something like PixelsToNativeRatio).

Error in validate.CleanGeometry on fractional coordinates

This is copied from tegola #338

---

The test case for this can be found here.

To run this:

• checkout the intersect_benchmark branch from my repo
• cd tegola/maths/validate
• Run: go test -v .

At this time the test is only printing results and isn't flagging the tests as failures.

Input polygon and extent:

// fractional coordinate
// if you multiply everything by 10 (remove fraction), then this test case
// works just fine.
{
	basic.Polygon { basic.Line {
			basic.Point {.1, 0},
			basic.Point {1, 0},
			basic.Point {0, 1},
	} },
	geom.NewExtent([2]float64{0, 0}, [2]float64{.5, .5}),
},

The winding order and POLYGON vs. MULTIPOLYGON aren't really relevant, the relevant part is (0.05, .5) point vs (0, 0.5) point.

Expected: POLYGON ((0.5 0,0.1 0,0.05 0.5,0.5 0.5))
Actual: MULTIPOLYGON (((0 0.5,0.1 0,0.5 0,0.5 0.5)))

`TestMakeValid/makevalidTestCases_#6_circle_one_right` is failing.

TestMakeValid/makevalidTestCases_#6_circle_one_right is failing.

However, the debugger visual seems to indicate that it should be passing:

Figure out what the issue is and get the test passing.

Slight decode behavior difference between wkt.Decode and geojson.UnmarshalJSON

I found a slight behavior diffrenece between 2 decoding methods.

• wkt.Decode removes a last overlapping point.
• geojson.UnmarshalJSON does not.

Is this known or an as-intended behavior?

-- EDIT
I forgot to write package's version.

uname -a

Linux 1b647680731a 5.10.60.1-microsoft-standard-WSL2 #1 SMP Wed Aug 25 23:20:18 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux

go version

go version go1.19 linux/amd64

geom

github.com/go-spatial/geom v0.0.0-20220918193402-3cd2f5a9a082

// https://go.dev/play/p/zX8V3DC9N2_I
// You can edit this code!
// Click here and start typing.
package main

import (
	"encoding/json"
	"fmt"

	"github.com/go-spatial/geom"
	"github.com/go-spatial/geom/encoding/geojson"
	"github.com/go-spatial/geom/encoding/wkt"
)

func main() {
	geojsonBin := []byte(`{
			"type" : "Polygon",
			"coordinates" : [[
					 [100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0]
			]]
		}`)
	wktText := "POLYGON ((100.0 0.0, 101.0 0.0, 101.0 1.0, 100.0 1.0, 100.0 0.0))"

	var jsonG geojson.Geometry
	json.Unmarshal(geojsonBin, &jsonG)

	fmt.Println(jsonG.Geometry.(geom.Polygon).LinearRings())
	// It does not remove the last overlapping point.
	// [[[100 0] [101 0] [101 1] [100 1] [100 0]]]

	wktG, _ := wkt.DecodeString(wktText)
	fmt.Println(wktG.(geom.Polygon).LinearRings())
	// It removes the last point.
	// [[[100 0] [101 0] [101 1] [100 1]]]
}

[geom] fix NaN comparisons

For ieee754 floats comparing NaN with NaN is always false:

https://play.golang.org/p/UbVEpc9zaSR

We can either use x != x to find them or use the math.IsNaN library function which does exactly that:

https://golang.org/src/math/bits.go?s=791:822#L24

Here's one place in the codebase that compares NaN with each other:

geom/geom.go

Line 507 in a43924e

return g[0] == nan && g[1] == nan

geom/point.go

Line 11 in a43924e

var nan = math.NaN()

typo in LineString.Verticies

The plural of vertex is vertices, however there might be dependents of this package. Should we change this?

https://www.merriam-webster.com/dictionary/vertex

Panic in the edge_triangle.go.

When cutting out an internal "ab … bc" ring, we get a panic.

2018/08/02 17:38:20 edge_triangle.go:345: oog rng: [[985 1485] [986 1484] [986 1483] [989 1483] [991 1482] [993 1483] [993 1485] [992 1487] [989 1489] [988 1490] [987 1489] [986 1487] [988 1485] [986 1487]]
2018/08/02 17:38:20 edge_triangle.go:346: org rng: [[985 1485] [986 1483] [989 1483] [991 1482] [993 1483] [993 1485] [992 1487] [989 1489] [988 1490] [987 1489]]
2018/08/02 17:38:20 edge_triangle.go:347: bubble type ab…bc: (    2)(    3) … (    2)(    3)
2018/08/02 17:38:20 edge_triangle.go:348: bubble type ab…bc: ([989 1483])([991 1482]) … ([989 1483])([991 1482])
2018/08/02 17:38:20 edge_triangle.go:349: cur rng: [[989 1483]]
2018/08/02 17:38:20 edge_triangle.go:350: sliver : [[991 1482] [993 1483] [993 1485] [992 1487] [989 1489] [988 1490] [987 1489] [985 1485] [986 1483]]
2018/08/02 17:38:20 edge_triangle.go:345: oog rng: [[985 1485] [986 1484] [986 1483] [989 1483] [991 1482] [993 1483] [993 1485] [992 1487] [989 1489] [988 1490] [987 1489] [986 1487] [988 1485] [986 1487]]
2018/08/02 17:38:20 edge_triangle.go:346: org rng: [[985 1485] [986 1484] [986 1483] [989 1483] [991 1482] [993 1483] [993 1485] [992 1487] [989 1489] [988 1490] [987 1489] [986 1487] [988 1485] [986 1487]]
2018/08/02 17:38:20 edge_triangle.go:347: bubble type ab…bc: (   10)(    3) … (   13)(    0)
2018/08/02 17:38:20 edge_triangle.go:348: bubble type ab…bc: ([987 1489])([989 1483]) … ([986 1483])([985 1485])
2018/08/02 17:38:20 edge_triangle.go:349: cur rng: [[989 1483]]
2018/08/02 17:38:20 edge_triangle.go:350: sliver : [[986 1487] [988 1485]]
2018/08/02 17:38:20 server.go:2923: http: panic serving [::1]:60302: index out of bounds[0 - 1], start 3 end 4

[encoding] WKT

In the latest version of the WKT code, I did not do a good job of code reviewing it.

After using it and reading the code base a bit more, I have a few issues I feel that should be fixed:

1. move isNil to cmp package.

2. move IsEmptyPoint to the cmp package

Both these functions are useful functions for more than this package.

1. Encoder struct should have two members:
   1.1. Strict bool -- this should be a bool that changes the encoder's behavior so that if Linestring, Multilinestrings, Multi points, Polygons, and Multipolygons mix empty points and non_empty points in their make up the encoder will error. The default should be false.
   1.1.1. If strict is false, then Linestring, Multilinestring, Multipoints, Polygons, and Multipolygons with mixed points, will drop the empty points.
   1.2. precision this should be the value that is used here. QGIS does not support the 'e' notation. So, I had to bump this up to 12. I think 5 is too small anyway.

Rename isPointless to something that makes more sense.
For example rename isPointlessPoint to lastIndexEmptyPointForPoints(mp [][2]float64)(idx int), where idx is -1 if no empty points were found.

Rename removePointless to removeEmptyPoints params for various functions.
Though that should not be passed in anymore, and instead just be the default behavior. With Strict reporting if there is mixed points being used as apporiate.

I know putc and puts are prettty standard c functions. but it's hard to understand what is going on. I like just character and string and the reciever is already enc (encode) it reads better: enc.character('c')

Polygon cut is appearing as an additional polygon

As you can see from this picture: (location of the issue: http://localhost:8080/?debug=true#6.62/23.78/-76.40)

The pointed out outline should be a cutout, and not solid. This indicates that the cutout is being treated as a seperate polygon

Provided are three wkt's of the region.

This is the full wkt
geom_query.csv.txt

This is the area with in the extent st_makeenvelope(-8.77619383836914e+06,2.4949046028808597e+06,-8.13045382350586e+06,3.140644617744142e+06,3857)

geom_query_got.csv.txt

One that is focused in on the issue:
geom_query_mini.csv.txt

[ci] CGO and No CGO being compared in coverage tests

We should either turn off one of them (No CGO has my vote), or figure out how to separate them in coverall api

To disable push of the cover information we just need to remove the --coveralls. To have different coverprofile, we just need to pass in --coverprofile=$name

The easiest thing might be to use the `--coverprofile=default${CGO_ENABLED}" flag to the coverall and that should have different cover profiles?

But this is a bit of a hack.

Grid: Unpredictable interactions between ToNative and FromNative

Expected

Grid.ToNative() should return the top-left point of the given tile:

grid, _ := slippy.NewGrid(3857)
tile := slippy.NewTile(5, 0, 0)
point, _ := grid.ToNative(tile)
fmt.Println(point) // [-2.003750834e+07 2.003750834e+07]

Grid.FromNative() should return the tile corresponding to a top-left point:

grid, _ := slippy.NewGrid(3857)
tile, _ := grid.FromNative(5, geom.Point{-2.003750834e7, 2.003750834e7})
fmt.Println(tile) // &{5 0 0}

FromNative() and ToNative() should be able to reverse each other predictably. That is, FromNative(ToNative(x)) == x.

Problem

The above does not hold for all coordinates. In the following example, 5,3,5 becomes 5,3,4:

grid, _ := NewGrid(3857)
tile1 := NewTile(5, 3, 5)
point, _ := grid.ToNative(tile1)
tile2, _ := grid.FromNative(5, point)
fmt.Println(tile2) // &{5 3 4}

x=3 and y=5 seems to be buggy at multiple zooms (try 6,3,5, 7,3,5, and 8,3,5)

Impact

Whatever is causing this problem seems to make slippy.RangeFamilyAt() unpredictable as well, as it relies on ToNative() and FromNative(). When updating Tegola to use the latest geom (go-spatial/tegola#952), I couldn't get RangeFamilyAt() to work and had to copy in the old version which does not make use of these functions.

Note: it's entirely possible this is related to me using an M1 Mac. I see slightly smaller floats generated on a regular basis, so maybe someone could try this on amd64?

[triangulation] Failed to insert point

Failed to insert point while triangulation.

https://gist.github.com/gdey/38eca8aea9ae85e5d2c7f8946d89f910

---

Possible additional test cases for triangulations:
https://gist.github.com/ear7h/b231f1bf4ff4b3c74555a0c2fb30208d

mvt: add README.md

Add README.md file to the mvt package.

for gpkg encoding library.

Currently the gpkg only support what is needed by us.

However, here is the sql from QGIS, there are some interesting triggers for the tile matrix:

			allSQLFROMQGIS = `

		CREATE TABLE IF NOT EXISTS gpkg_spatial_ref_sys  (srs_name TEXT NOT NULL,srs_id INTEGER NOT NULL PRIMARY KEY,organization TEXT NOT NULL,organization_coordsys_id INTEGER NOT NULL,definition  TEXT NOT NULL,description TEXT);
		CREATE TABLE IF NOT EXISTS gpkg_contents  (table_name TEXT NOT NULL PRIMARY KEY,data_type TEXT NOT NULL,identifier TEXT UNIQUE,description TEXT DEFAULT '',last_change DATETIME NOT NULL DEFAULT (strftime('%Y-%m-%dT%H:%M:%fZ','now')),min_x DOUBLE, min_y DOUBLE,max_x DOUBLE, max_y DOUBLE,srs_id INTEGER,CONSTRAINT fk_gc_r_srs_id FOREIGN KEY (srs_id) REFERENCES gpkg_spatial_ref_sys(srs_id));
		CREATE TABLE IF NOT EXISTS gpkg_ogr_contents (table_name TEXT NOT NULL PRIMARY KEY,feature_count INTEGER DEFAULT NULL);
		CREATE TABLE IF NOT EXISTS gpkg_geometry_columns (table_name TEXT NOT NULL,column_name TEXT NOT NULL,geometry_type_name TEXT NOT NULL,srs_id INTEGER NOT NULL,z TINYINT NOT NULL,m TINYINT NOT NULL,CONSTRAINT pk_geom_cols PRIMARY KEY (table_name, column_name),CONSTRAINT uk_gc_table_name UNIQUE (table_name),CONSTRAINT fk_gc_tn FOREIGN KEY (table_name) REFERENCES gpkg_contents(table_name),CONSTRAINT fk_gc_srs FOREIGN KEY (srs_id) REFERENCES gpkg_spatial_ref_sys (srs_id));
		CREATE TABLE IF NOT EXISTS gpkg_tile_matrix_set (table_name TEXT NOT NULL PRIMARY KEY,srs_id INTEGER NOT NULL,min_x DOUBLE NOT NULL,min_y DOUBLE NOT NULL,max_x DOUBLE NOT NULL,max_y DOUBLE NOT NULL,CONSTRAINT fk_gtms_table_name FOREIGN KEY (table_name) REFERENCES gpkg_contents(table_name),CONSTRAINT fk_gtms_srs FOREIGN KEY (srs_id) REFERENCES gpkg_spatial_ref_sys (srs_id));
		CREATE TABLE IF NOT EXISTS gpkg_tile_matrix (table_name TEXT NOT NULL,zoom_level INTEGER NOT NULL,matrix_width INTEGER NOT NULL,matrix_height INTEGER NOT NULL,tile_width INTEGER NOT NULL,tile_height INTEGER NOT NULL,pixel_x_size DOUBLE NOT NULL,pixel_y_size DOUBLE NOT NULL,CONSTRAINT pk_ttm PRIMARY KEY (table_name, zoom_level),CONSTRAINT fk_tmm_table_name FOREIGN KEY (table_name) REFERENCES gpkg_contents(table_name));
		CREATE TRIGGER IF NOT EXISTS 'gpkg_tile_matrix_zoom_level_insert' BEFORE INSERT ON 'gpkg_tile_matrix' FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'insert on table ''gpkg_tile_matrix'' violates constraint: zoom_level cannot be less than 0') WHERE (NEW.zoom_level < 0); END;
		CREATE TRIGGER IF NOT EXISTS 'gpkg_tile_matrix_zoom_level_update' BEFORE UPDATE of zoom_level ON 'gpkg_tile_matrix' FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'update on table ''gpkg_tile_matrix'' violates constraint: zoom_level cannot be less than 0') WHERE (NEW.zoom_level < 0); END;
		CREATE TRIGGER IF NOT EXISTS 'gpkg_tile_matrix_matrix_width_insert' BEFORE INSERT ON 'gpkg_tile_matrix' FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'insert on table ''gpkg_tile_matrix'' violates constraint: matrix_width cannot be less than 1') WHERE (NEW.matrix_width < 1); END;
		CREATE TRIGGER IF NOT EXISTS 'gpkg_tile_matrix_matrix_width_update' BEFORE UPDATE OF matrix_width ON 'gpkg_tile_matrix' FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'update on table ''gpkg_tile_matrix'' violates constraint: matrix_width cannot be less than 1') WHERE (NEW.matrix_width < 1); END;
		CREATE TRIGGER IF NOT EXISTS 'gpkg_tile_matrix_matrix_height_insert' BEFORE INSERT ON 'gpkg_tile_matrix' FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'insert on table ''gpkg_tile_matrix'' violates constraint: matrix_height cannot be less than 1') WHERE (NEW.matrix_height < 1); END;
		CREATE TRIGGER IF NOT EXISTS 'gpkg_tile_matrix_matrix_height_update' BEFORE UPDATE OF matrix_height ON 'gpkg_tile_matrix' FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'update on table ''gpkg_tile_matrix'' violates constraint: matrix_height cannot be less than 1') WHERE (NEW.matrix_height < 1); END;
		CREATE TRIGGER IF NOT EXISTS 'gpkg_tile_matrix_pixel_x_size_insert' BEFORE INSERT ON 'gpkg_tile_matrix' FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'insert on table ''gpkg_tile_matrix'' violates constraint: pixel_x_size must be greater than 0') WHERE NOT (NEW.pixel_x_size > 0); END;
		CREATE TRIGGER IF NOT EXISTS 'gpkg_tile_matrix_pixel_x_size_update' BEFORE UPDATE OF pixel_x_size ON 'gpkg_tile_matrix' FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'update on table ''gpkg_tile_matrix'' violates constraint: pixel_x_size must be greater than 0') WHERE NOT (NEW.pixel_x_size > 0); END;
		CREATE TRIGGER IF NOT EXISTS 'gpkg_tile_matrix_pixel_y_size_insert' BEFORE INSERT ON 'gpkg_tile_matrix' FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'insert on table ''gpkg_tile_matrix'' violates constraint: pixel_y_size must be greater than 0') WHERE NOT (NEW.pixel_y_size > 0); END;
		CREATE TRIGGER IF NOT EXISTS 'gpkg_tile_matrix_pixel_y_size_update' BEFORE UPDATE OF pixel_y_size ON 'gpkg_tile_matrix' FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'update on table ''gpkg_tile_matrix'' violates constraint: pixel_y_size must be greater than 0') WHERE NOT (NEW.pixel_y_size > 0); END;

			`

cmp: epsilon not suitble for large numbers

The issue

The core of the cmp package is the following lines

geom/cmp/cmp.go

Line 62 in c4fa59f

return math.Abs(f1-f2) < tolerance

Where tolerance is generally the package constant TOLERANCE. This approach eases some of the problems of floats, but it breaks if numbers are too large or too small. For example, I am currently working with points in the Natural Earth dataset, the the following comparison fails despite being accurate to 10 significant digits.

cmp.Float(
    -7.784854276e+06,
    -7.78485427595008e+06) // false

proposed solution

The following page provides more details and a solution:

https://floating-point-gui.de/errors/comparison/

There is an alternative to heaping conceptual complexity onto such an apparently simple
task: instead of comparing a and b as real numbers, we can think about them as discrete
steps and define the error margin as the maximum number of possible floating-point values
between the two values.

This is conceptually very clear and easy and has the advantage of implicitly scaling the
relative error margin with the magnitude of the values. Technically, it’s a bit more complex,
but not as much as you might think, because IEEE 754 floats are designed to maintain their
order when their bit patterns are interpreted as integers.

In code, this would be something like:

        // comparing against zero needs extra complexity bc theres
        // the just as many bitpatterns between [0, 1] as [1, inf)
        // and the count method breaks down. So, we turn it into a comparison
        // in the [1, inf) range.
        if f1 == 0 {
                f1 = 1
                f2 = math.Abs(f2) + 1
        }
        if f2 == 0 {
                f1 = math.Abs(f1) + 1
                f2 = 1
        }

        i1 := *(*int64)(unsafe.Pointer(&f1))
        i2 := *(*int64)(unsafe.Pointer(&f2))
        d := i2 - i1
        if d < 0 {
                return d > -tolerance
        } else {
                return d < tolerance
        }

The new tolerance constant could be calculated as the number of steps between 1.000001 and 1.0, which is 4503599627 (see this play). This would make the comparison accurate to 5 decimal places.

[makevalid] [constrainted triangulation] panic in (*QuadEdge).Sym

Testing the new simplification function, I am running into a nil pointer dereference.

panic: runtime error: invalid memory address or nil pointer dereference [recovered]
	panic: runtime error: invalid memory address or nil pointer dereference
[signal SIGSEGV: segmentation violation code=0x1 addr=0x0 pc=0x467a656]

goroutine 407 [running]:
github.com/go-spatial/tegola/atlas.Map.Encode.func1.1.1(0xc0002e1aa0, 0xc000c67e00, 0xc0070f7358)
	/Users/julio/lib/go/src/github.com/go-spatial/tegola/atlas/map.go:215 +0x505
panic(0x49c0580, 0x527c920)
	/usr/local/go/src/runtime/panic.go:679 +0x1b2
github.com/go-spatial/geom/planar/triangulate/quadedge.(*QuadEdge).Sym(...)
	/Users/julio/lib/go/src/github.com/go-spatial/geom/planar/triangulate/quadedge/quadedge.go:251
github.com/go-spatial/geom/planar/triangulate/constraineddelaunay.(*Triangulator).deleteEdge(0xc0070f6f58, 0xc0009644e0, 0x41716c4ca9fb7e18, 0x415ffcf349193de5)
	/Users/julio/lib/go/src/github.com/go-spatial/geom/planar/triangulate/constraineddelaunay/triangulator.go:162 +0xa6
github.com/go-spatial/geom/planar/triangulate/constraineddelaunay.(*Triangulator).insertEdgeCDT(0xc0070f6f58, 0xc0070f6a88, 0x0, 0x0, 0x0, 0x0)
	/Users/julio/lib/go/src/github.com/go-spatial/geom/planar/triangulate/constraineddelaunay/triangulator.go:757 +0x1678
github.com/go-spatial/geom/planar/triangulate/constraineddelaunay.(*Triangulator).insertConstraints(0xc0070f6f58, 0x49d6e80, 0xc0014c3320, 0x0, 0x0, 0x0, 0x42d7c700ad4cb1c5)
	/Users/julio/lib/go/src/github.com/go-spatial/geom/planar/triangulate/constraineddelaunay/triangulator.go:460 +0x5a8
github.com/go-spatial/geom/planar/triangulate/constraineddelaunay.(*Triangulator).InsertGeometries(0xc0070f6f58, 0xc0070f6f48, 0x1, 0x1, 0x0, 0x0, 0x0, 0xc0070f6ee0, 0x466d202)
	/Users/julio/lib/go/src/github.com/go-spatial/geom/planar/triangulate/constraineddelaunay/triangulator.go:373 +0x13b
github.com/go-spatial/geom/planar/makevalid.TriangulateGeometry(0x4c38980, 0xc0008cd350, 0x49d6e80, 0xc0014c3320, 0x54, 0x5c, 0x56, 0xc0070f7010, 0x409f999)
	/Users/julio/lib/go/src/github.com/go-spatial/geom/planar/makevalid/triangulate.go:19 +0xb3
github.com/go-spatial/geom/planar/makevalid.InsideTrianglesForGeometry(0x4c38980, 0xc0008cd350, 0x49d6e80, 0xc0014c3320, 0x4c31c00, 0xc00411fd20, 0x0, 0x0, 0x0, 0xc000099800, ...)
	/Users/julio/lib/go/src/github.com/go-spatial/geom/planar/makevalid/triangulate.go:57 +0x6a
github.com/go-spatial/geom/planar/makevalid.(*Makevalid).makevalidPolygon(0xc0070f7388, 0x4c38980, 0xc0008cd350, 0xc000798020, 0xc0014c3320, 0x1, 0x0, 0xc0070f7218)
	/Users/julio/lib/go/src/github.com/go-spatial/geom/planar/makevalid/makevalid.go:194 +0x11c
github.com/go-spatial/geom/planar/makevalid.(*Makevalid).Makevalid(0xc0070f7388, 0x4c38980, 0xc0008cd350, 0x49cd360, 0xc00683a2a0, 0xc000798020, 0x0, 0xc00683a2a0, 0x0, 0x0, ...)
	/Users/julio/lib/go/src/github.com/go-spatial/geom/planar/makevalid/makevalid.go:235 +0x57c
github.com/go-spatial/tegola/atlas.Map.Encode.func1.1(0xc0015a0000, 0x0, 0x0)
	/Users/julio/lib/go/src/github.com/go-spatial/tegola/atlas/map.go:254 +0x384
github.com/go-spatial/tegola/provider/postgis.Provider.TileFeatures(0xc0000d1570, 0x9, 0x1538, 0xc0000d15a0, 0xd, 0xc0000d15d0, 0x6, 0x0, 0x0, 0x0, ...)
	/Users/julio/lib/go/src/github.com/go-spatial/tegola/provider/postgis/postgis.go:601 +0xa41
github.com/go-spatial/tegola/atlas.Map.Encode.func1(0xc000cbe450, 0xc000c67e00, 0xc000255d90, 0x2, 0xc00027ede0, 0x10, 0x520d660, 0xc053119efffcb851, 0x404393a5a4fa85e7, 0x4014000000000000, ...)
	/Users/julio/lib/go/src/github.com/go-spatial/tegola/atlas/map.go:166 +0x272
created by github.com/go-spatial/tegola/atlas.Map.Encode
	/Users/julio/lib/go/src/github.com/go-spatial/tegola/atlas/map.go:154 +0x208

slippy.NewTile should check tile if x and y are valid for given z

geom/slippy/tile.go

Lines 15 to 21 in 2760310

	func NewTile(z, x, y uint) *Tile {
	return &Tile{
	Z: z,
	X: x,
	Y: y,
	}
	}

The NewTile function simply fills in the fields but it should do input validation, for example the tile {Z:0, X: 99, Y:99} does not exist.

I can see two solutions, changing the function signature to func NewTile(z, x, y uint) (*Tile, error) or returning nil from the current function if the tile does not exist.

New Signature

pros:

• more explicit error state
• error state is a compile error
  • tile := slippy.NewTile(...); tile.Dereference is caught by the compiler
• the returned error can be created in geom

cons:

• breaking change

Return nil

pros:

• non-breaking
• error state would give a runtime error (better than going un-noticed)
  • tile := slippy.NewTile(...); tile.Dereference is a nil pointer dereference at runtime

cons:

• error state is would give a runtime error (not as helpful as a compiler error)