diff --git a/README.md b/README.md index f31ed18..97e986d 100644 --- a/README.md +++ b/README.md @@ -7,69 +7,52 @@ Card](https://goreportcard.com/badge/PDOK/texel)](https://goreportcard.com/repor [](https://hub.docker.com/r/pdok/texel) -Preprocesses a [GeoPackage](https://www.geopackage.org/) -for creating _valid_ (more info below) vector tiles from it. -(And in the future probably will also create the actual vector tiles.) - -The (MULTI)POLYGON geometries in a geopackage are snapped to the grid cq matrix -matching the one of the intended (internal pixel grid of the) vector tiles. -While ensuring that no (new) intersections arise from this, by adding points -to lines that skim along other vertices. -An additional benefit is that the vector tile generator needs to process fewer -data (snapping combines some points into one) and does not need any additional -processing aside from clipping and converting to the internal tile pixel grid. +Processes a [GeoPackage](https://www.geopackage.org/) in order to produce Mapbox vector tiles. -## "Valid" geometries - -Our definition of _valid_ is in the context of vector tiles and is more loose -than the [OGC rules](https://en.wikipedia.org/wiki/Simple_Features). -Mainly in comparison to the OGC rules, we allow: - -* Overlap (not intersection) -* Linear rings or line strings with less than 3 vertices - (effectively becoming lines or points) - -With the drawing/rendering of vector tiles in mind, these invalidities do not -cause a problem. On the contrary, they introduce possibilities. -E.g. when zooming out, you would like to still see a lake as a point and a -river as a line. (The tool that this codebase is forked from, `sieve` would -filter those out.) Also when zooming out on a lake with an attached river, -you would like to keep seeing that river as a line, but still as part of the -original geometry/feature including all attributes and styling. +This tool was originally created to snap MULTI(POLYGONS) to the grid cq matrix +corresponding to the vector tiles, while preserving a weak notion of _valid_ +(more info below). Subsequently, the tool can be used in two ways: +- `texel snap`, which snaps the polygons. +- `texel snap -enc` with `texel mvt`. The first command snaps and encodes + geometries according to the standard. The second assembles the tiles. -## Examples - -Not adding extra points could create intersections: - - - -Close thin polygons can turn into lines: - -  +## Usage -## Notes +Getting help: -- It will take a Geopackage and writes a new Geopackage with all the - (MULTI)POLYGON tables preprocessed. - - All other spatial tables are 'untouched' and copied as-is. - - Other non-spatial tables are not copied to the new geopackage. -- :warning: Spatialite lib is mandatory for running this application. This lib is needed for - creating the RTree triggers on the spatial tables for updating/maintaining the - RTree. +```sh +./texel --help +./texel snap --help +./texel mvt --help +``` -## Usage +Snap and encode geopackage ```sh -./texel -s=[source GPKG] -t=[target GPKG] \ +./texel snap \ + -s=[source GPKG] -t=[target GPKG] \ -tms=[tile matrix set for filtering] -z=[tile matrix ids] \ -p=[pagesize for writing to target GPKG] -o=[overwrite target GPKG] \ - -pl=[keep points and lines] + -pl=[keep points and lines] \ + -iog=[ignore outside grid] --rwo=[reverse winding order] \ + -enc=[encode geometries] \ + -buf=[buffer size for tile detection] \ + -lt=[line trace method for tile detection] \ + -c=[current clip] +``` -./texel --help +Make tiles from snapped output + +```sh +./texel mvt \ + -c=[config toml file] \ + -z=[tile matrix id] ``` ### Docker +Example Docker command for running texel without encoding geometries. + ```docker docker run \ --name texel \ @@ -77,6 +60,7 @@ docker run \ -u $(id -u):$(id -g) \ -v `pwd`/example:/example \ pdok/texel \ + snap -s=./example/example.gpkg \ -t=./example/example-processed.gpkg \ -tms="NetherlandsRDNewQuad" \ @@ -100,9 +84,152 @@ go build . docker build -t pdok/texel . ``` +:warning: Spatialite lib is mandatory for running this application. This lib is +needed for creating the RTree triggers on the spatial tables for +updating/maintaining the RTree. + +## Commands + +### Snapping + +`texel snap` takes a geopackage, snaps all (MULTI)POLYGON geometries in it, and +writes to a new geopackage. If encoding is enabled, for each table it also +creates `{TABLENAME}_encoded` tables. A row in such a table is a geometry +encoded for a particular tile. These are required input for `texel mvt` +(explained below). + +- All other spatial tables are 'untouched' and copied as-is. +- Other non-spatial tables are not copied to the new geopackage. + +Input parameters: +- source (`-s`): source Geopackage. +- target (`-t`): target Geopackage. If specified as `target.gpkg`, output is + written to `target_{zoomlevel}.gpkg` for each zoomlevel. +- tile matrix set (`-tms`): tile matrix set name. Supported tile matrix sets can + be found in the `tms20/tilematrixsets` folder. +- tile matrix ids (`-z`): JSON array of integers of zoomlevels to target. + +Snap behaviour parameters: +- keep points and lines: (`-pl`): if (part of) a polygon degenerates to points + and lines during snapping remove those parts. +- ignore outside grid (`-iog`): skip geometries outside of the grid +- reverse winding order (`-rwo`): reverse winding order of geometries. This is + used to comply with the vectortile standard. + +Encoding behaviour parameters: +- encode (`-enc`): enable encoding. If set to false, do not create the + `{TABLENAME}_encoded` tables. +- line trace (`-lt`): there are two algorithms to determine on which tiles a + geometry lies. When `-lt=false`, use a bounding box. When `-lt=true`, trace + border of features. This is more precise, but more computationally expensive. +- buffer (`-buf`): buffer size around each tile. Geometries within this buffer + are encoded for this tile. +- clip (`-c`): a higher-level tile written as `[z, x, y]`, for zoomlevel and + coordinates. It is required that this `z` not larger than the + `tile-matrix-ids` above. Only encode geometries lying within this + higher-level tile. + +The `clip` parameter might require more explanation. One possible workflow is +divide a workload of generating all tiles into several "clips", for example for +parallel processing. A clip is such a tile at a higher level. In many +instances, for example when the buffer is strictly positive, the standard +workflow of `texel snap` + `texel mvt` generates tiles neighbouring the clip. +This option prevents that. + +## Tile generation + +The output of `texel snap -enc` is input for `texel mvt`, which generates +vectortiles. There are only a few parameters + +- config (`-c`): a TOML configuration file that is detailed below. +- tile matrix id (`-z`): the zoomlevel that is being targeted. + +Almost all configuration is determined in this TOML file. The configuration is +similar to configurations for +[trex](https://github.com/t-rex-tileserver/t-rex), however, valid +configurations for `trex` are not valid for `texel` in a way detailed below. +The following configuration is an example. + +```toml +[tileset] +name = "NetherlandsRDNewQuad" + +[[datasource]] +name = "brt" +path = "texeled/brt_12.gpkg" + +[[tileset.layer]] +name = "wegdeel" +datasource = "brt" +table_name = "wegdeel_12" +minzoom = 12 +maxzoom = 12 + +[[tileset.layer]] +name = "waterdeel" +datasource = "brt" +table_name = "waterdeel_11_12" +minzoom = 11 +maxzoom = 12 + +[cache.file] +base = "/data" + +[cache.azure] +connection_string = "{{env.AZURE_CONNECTION_STRING}}" +container = "{{env.AZURE_CONTAINER}}" +prefix_key = "{{env.AZURE_KEY_PREFIX}}" +``` + +Notes: +- A `[[datasource]]` is a geopackage that is results from `texel snap -enc`. +- A `[[tileset.layer]]` represents both a table in a geopackage and a layer in + the resulting vectortile. The geopackage and table is determined by + `datasource`, which refers to the name of a datasource and `table_name` + (which should NOT inclde an `_encoded` suffix). + + A layer is included if the parameter `-z` specified above lies between + `minzoom` and `maxzoom` (inclusive). +- If `[cache.file]` is specified, tiles are written with that path as base. +- If `[cache.azure]` is specified, write to the an azure blobstore. +- At least one of `[cache.file]` or `[cache.azure]` needs to be specified. +- At the destination, the tiles are written as `{NAME}/{z}/{x}/{y}.pbf`. Here + `name` is the name of the tileset. When writing to the azure blobstore, the + `prefix_key` is also prepended. + +## "Valid" geometries + +`texel snap` produces "valid" geometries. + +Our definition of _valid_ is in the context of vector tiles and is more loose +than the [OGC rules](https://en.wikipedia.org/wiki/Simple_Features). +Mainly in comparison to the OGC rules, we allow: + +* Overlap (not intersection) +* Linear rings or line strings with less than 3 vertices + (effectively becoming lines or points) +* Gaps + +With the drawing/rendering of vector tiles in mind, these invalidities do not +cause a problem. On the contrary, they introduce possibilities. +E.g. when zooming out, you would like to still see a lake as a point and a +river as a line. (The tool that this codebase is forked from, `sieve` would +filter those out.) Also when zooming out on a lake with an attached river, +you would like to keep seeing that river as a line, but still as part of the +original geometry/feature including all attributes and styling. + +### Examples +Not adding extra points could create intersections: + + + +Close thin polygons can turn into lines: + +  ## References * [sieve](https://github.com/pdok/sieve) is the predecessor of this tool and much of its processing codebase is reused as boilerplate. + diff --git a/config/mvtconfig.go b/config/mvtconfig.go new file mode 100644 index 0000000..2b27ed4 --- /dev/null +++ b/config/mvtconfig.go @@ -0,0 +1,101 @@ +package config + +import ( + "fmt" + "os" + "regexp" + + "github.com/BurntSushi/toml" +) + +type DataSource struct { + Name string `toml:"name"` + Path string `toml:"path"` +} + +type LayerConfig struct { + Name string `toml:"name"` + MinZoom uint `toml:"minzoom"` + MaxZoom uint `toml:"maxzoom"` + DataSource string `toml:"datasource"` + TableName string `toml:"table_name"` +} + +// Tileset mirrors a `[[Tileset]]` array-of-tables entry, including its +// nested `[[Tileset.layer]]` entries. +type Tileset struct { + Name string `toml:"name"` + Layer []LayerConfig `toml:"layer"` +} + +type AzureConfig struct { + ConnectionString string `toml:"connection_string"` + Container string `toml:"container"` + PrefixKey string `toml:"prefix_key"` +} + +type FileConfig struct { + Base string `toml:"base"` +} + +// In case S3 will be needed in the future, these are the old `t_rex` configs. +// type S3Config struct { +// Endpoint string `toml:"endpoint"` +// Bucket string `toml:"bucket"` +// AccessKey string `toml:"access_key"` +// SecretKey string `toml:"secret_key"` +// Region string `toml:"region"` +// KeyPrefix string `toml:"key_prefix"` +// } + +type CacheConfig struct { + File *FileConfig `toml:"file"` + Azure *AzureConfig `toml:"azure"` +} + +// TomlConfig mirrors the top-level structure of an mvt config toml file, +// e.g. example/NetherlandsRDNewQuad.toml. +type TomlConfig struct { + DataSource []DataSource `toml:"datasource"` + Tileset *Tileset `toml:"tileset"` + Cache CacheConfig `toml:"cache"` +} + +// Replace {{env.NAME}} by value of environment variable NAME. Inserted +// verbatim. Unset variables result in an error. +func substituteEnv(s string) (string, error) { + envPlaceholder := regexp.MustCompile( + `\{\{\s*env\.([A-Za-z_][A-Za-z0-9_]*)\s*\}\}`) + + var missing []string + out := envPlaceholder.ReplaceAllStringFunc(s, func(m string) string { + name := envPlaceholder.FindStringSubmatch(m)[1] + val, ok := os.LookupEnv(name) + if !ok { + missing = append(missing, name) + } + return val + }) + if len(missing) > 0 { + return "", fmt.Errorf("unset environment variables: %v", missing) + } + return out, nil +} + +// ParseMVTConfig reads an mvt config toml file, substitutes {{env.NAME}} +// placeholders and returns the resulting TomlConfig. +func ParseMVTConfig(path string) (TomlConfig, error) { + var cfg TomlConfig + raw, err := os.ReadFile(path) + if err != nil { + return TomlConfig{}, fmt.Errorf("reading mvt config %q: %w", path, err) + } + text, err := substituteEnv(string(raw)) + if err != nil { + return TomlConfig{}, fmt.Errorf("mvt config %q: %w", path, err) + } + if _, err := toml.Decode(text, &cfg); err != nil { + return TomlConfig{}, fmt.Errorf("decoding mvt config %q: %w", path, err) + } + return cfg, nil +} diff --git a/config/mvtconfig_test.go b/config/mvtconfig_test.go new file mode 100644 index 0000000..e677eba --- /dev/null +++ b/config/mvtconfig_test.go @@ -0,0 +1,40 @@ +package config + +import ( + "testing" + + "github.com/stretchr/testify/require" +) + +func TestSubstituteEnv(t *testing.T) { + t.Setenv("TEXEL_TEST_PATH", "/data/a.gpkg") + t.Setenv("TEXEL_TEST_QUOTED", `\"TEST\"`) + + tests := []struct { + name string + in string + want string + }{ + {"no placeholder", `path = "x"`, `path = "x"`}, + {"simple", `path = "{{env.TEXEL_TEST_PATH}}"`, `path = "/data/a.gpkg"`}, + {"spaces", `path = "{{ env.TEXEL_TEST_PATH }}"`, `path = "/data/a.gpkg"`}, + { + "multiple", + `{{env.TEXEL_TEST_PATH}}:{{env.TEXEL_TEST_PATH}}`, + `/data/a.gpkg:/data/a.gpkg`, + }, + {"escaped quotes", `name = "{{env.TEXEL_TEST_QUOTED}}"`, `name = "\"TEST\""`}, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + got, err := substituteEnv(tt.in) + require.NoError(t, err) + require.Equal(t, tt.want, got) + }) + } +} + +func TestSubstituteEnv_Unset(t *testing.T) { + _, err := substituteEnv(`a = "{{env.TEXEL_TEST_DEFINITELY_UNSET}}"`) + require.ErrorContains(t, err, "TEXEL_TEST_DEFINITELY_UNSET") +} diff --git a/geomhelp/geomhelp.go b/geomhelp/geomhelp.go index 758a059..21a0930 100644 --- a/geomhelp/geomhelp.go +++ b/geomhelp/geomhelp.go @@ -1,6 +1,7 @@ package geomhelp import ( + "fmt" "math" "strings" @@ -146,7 +147,7 @@ func WktMustEncode(g geom.Geometry, maxLen uint) (s string) { return s } -func WktMustEncodeSlice(geoms []geom.Polygon, maxLen uint) string { +func WktMustEncodeSlice(geoms []geom.Geometry, maxLen uint) string { s := "" var builder strings.Builder for i := range geoms { @@ -162,3 +163,171 @@ func wktMustEncodeTruncated(geom geom.Geometry, width uint) string { } return truncate.StringWithTail(wkt.MustEncode(geom), width, "...") } + +// GeometrySliceToGeom converts a slice of points, linestrings or polygons +// into the matching (multi)geometry. Mixed geometry types are not allowed. +func GeometrySliceToGeom(geometries []geom.Geometry) geom.Geometry { + if len(geometries) == 0 { + return nil + } + if len(geometries) == 1 { + return geometries[0] + } + + switch geometries[0].(type) { + case geom.Point: + multiPoint := make(geom.MultiPoint, len(geometries)) + for i, g := range geometries { + p, ok := g.(geom.Point) + if !ok { + panic("Mixed geometry types are not allowed") + } + multiPoint[i] = p + } + return multiPoint + case geom.LineString: + multiLineString := make(geom.MultiLineString, len(geometries)) + for i, g := range geometries { + l, ok := g.(geom.LineString) + if !ok { + panic("Mixed geometry types are not allowed") + } + multiLineString[i] = l + } + return multiLineString + case geom.Polygon: + multiPolygon := make(geom.MultiPolygon, len(geometries)) + for i, g := range geometries { + p, ok := g.(geom.Polygon) + if !ok { + panic("Mixed geometry types are not allowed") + } + multiPolygon[i] = p + } + return multiPolygon + default: + panic("Unknown or unsupported geometry type") + } +} + +func PolygonSliceToGeom(polygons []geom.Polygon) geom.Geometry { + if len(polygons) == 0 { + panic("Multipolygon with zero polygons encountered") + } + if len(polygons) == 1 { + return polygons[0] + } + multipolygon := make(geom.MultiPolygon, len(polygons)) + for i, p := range polygons { + multipolygon[i] = p + } + return multipolygon +} + +// Convert geometry to slice of points. Does not accept multi-geometries +func PointSlice(g geom.Geometry) []geom.Point { + switch g := g.(type) { + case geom.Point: + return []geom.Point{g} + case geom.Polygon: + rings := g.LinearRings() + pointsCount := 0 + for _, ring := range rings { + pointsCount += len(ring) + } + points := make([]geom.Point, 0, pointsCount) + for _, ring := range rings { + for _, point := range ring { + points = append(points, point) + } + } + return points + case geom.LineString: + vertices := g.Vertices() + points := make([]geom.Point, 0, len(vertices)) + for _, vert := range vertices { + points = append(points, vert) + } + return points + default: + err := fmt.Errorf("PointSlice called incompatible geometry type: %s", g) + panic(err) + } +} + +// Essentially merge two lists +func mergeFamily[S any, E any, M ~[]E]( + g, h geom.Geometry, wrap func(S) M, errMsg string, +) geom.Geometry { + switch g := g.(type) { + case M: + switch h := h.(type) { + case M: + return append(g, h...) + case S: + return append(g, wrap(h)...) + default: + panic(errMsg) + } + case S: + switch h := h.(type) { + case M: + return append(wrap(g), h...) + case S: + return append(wrap(g), wrap(h)...) + default: + panic(errMsg) + } + default: + panic(errMsg) + } +} + +func MergeGeometries(g, h geom.Geometry) geom.Geometry { + if g == nil { + return h + } + if h == nil { + return g + } + + switch g.(type) { + case geom.Point, geom.MultiPoint: + return mergeFamily(g, h, func(p geom.Point) geom.MultiPoint { return geom.MultiPoint{p} }, + "Trying to merge point with non-point geometry") + case geom.LineString, geom.MultiLineString: + return mergeFamily(g, h, func(l geom.LineString) geom.MultiLineString { return geom.MultiLineString{l} }, + "Trying to merge linestring with non-linestring geometry") + case geom.Polygon, geom.MultiPolygon: + return mergeFamily(g, h, func(p geom.Polygon) geom.MultiPolygon { return geom.MultiPolygon{p} }, + "Trying to merge a polygon with a non-polygon geometry.") + default: + panic("Trying to merge unknown geometry types.") + } +} + +// A rather unfortunate function that needs to exist to deal with golang interfaces. +func MultiGeometryToSlice(multiGeom geom.Geometry) []geom.Geometry { + switch geometry := multiGeom.(type) { + case geom.MultiPolygon: + geoms := make([]geom.Geometry, 0, len(geometry)) + for _, polygon := range geometry { + geoms = append(geoms, geom.Polygon(polygon)) + } + return geoms + case geom.MultiLineString: + geoms := make([]geom.Geometry, 0, len(geometry)) + for _, line := range geometry { + geoms = append(geoms, geom.LineString(line)) + } + return geoms + case geom.MultiPoint: + geoms := make([]geom.Geometry, 0, len(geometry)) + for _, point := range geometry { + geoms = append(geoms, geom.Point(point)) + } + return geoms + default: + panic("Cannot convert to slice: not a multigeometry.") + } +} diff --git a/geomhelp/geomhelp_test.go b/geomhelp/geomhelp_test.go new file mode 100644 index 0000000..4085008 --- /dev/null +++ b/geomhelp/geomhelp_test.go @@ -0,0 +1,265 @@ +package geomhelp + +import ( + "testing" + + "github.com/go-spatial/geom" + "github.com/stretchr/testify/assert" +) + +func TestGeomhelp_PolygonSliceToGeom(t *testing.T) { + tests := []struct { + name string + polygonSlice []geom.Polygon + want geom.Geometry + wantPanic bool + }{ + { + name: "Empty polygon slice", + polygonSlice: []geom.Polygon{}, + wantPanic: true, + }, + { + name: "Single polygon", + polygonSlice: []geom.Polygon{{{{0.0, 0.0}, {1.0, 0.0}, {1.0, 1.0}}}}, + want: geom.Polygon{{{0.0, 0.0}, {1.0, 0.0}, {1.0, 1.0}}}, + }, + { + name: "Two polygons", + polygonSlice: []geom.Polygon{ + {{{0.0, 0.0}, {1.0, 0.0}, {1.0, 1.0}}}, + {{{10.0, 10.0}, {10.0, 11.0}, {11.0, 11.0}}}, + }, + want: geom.MultiPolygon{ + {{{0.0, 0.0}, {1.0, 0.0}, {1.0, 1.0}}}, + {{{10.0, 10.0}, {10.0, 11.0}, {11.0, 11.0}}}, + }, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + if tt.wantPanic { + assert.Panics(t, func() { PolygonSliceToGeom(tt.polygonSlice) }) + return + } + got := PolygonSliceToGeom(tt.polygonSlice) + assert.Equal(t, tt.want, got) + }) + } +} + +func TestMergeGeometries(t *testing.T) { + p1 := geom.Point{1, 2} + p2 := geom.Point{3, 4} + mp1 := geom.MultiPoint{{5, 6}, {7, 8}} + mp2 := geom.MultiPoint{{9, 9}} + + l1 := geom.LineString{{0, 0}, {1, 1}} + l2 := geom.LineString{{2, 2}, {3, 3}} + ml1 := geom.MultiLineString{{{4, 4}, {5, 5}}} + ml2 := geom.MultiLineString{{{6, 6}, {7, 7}}} + + poly1 := geom.Polygon{{{0, 0}, {1, 0}, {1, 1}}} + poly2 := geom.Polygon{{{2, 2}, {3, 2}, {3, 3}}} + mpoly1 := geom.MultiPolygon{{{{4, 4}, {5, 4}, {5, 5}}}} + mpoly2 := geom.MultiPolygon{{{{6, 6}, {7, 6}, {7, 7}}}} + + tests := []struct { + name string + g, h geom.Geometry + want geom.Geometry + }{ + {name: "nil g", g: nil, h: p1, want: p1}, + {name: "nil h", g: p1, h: nil, want: p1}, + + {name: "point + point", g: p1, h: p2, want: geom.MultiPoint{p1, p2}}, + {name: "point + multipoint", g: p1, h: mp1, want: geom.MultiPoint{{1, 2}, {5, 6}, {7, 8}}}, + {name: "multipoint + point", g: mp1, h: p1, want: geom.MultiPoint{{5, 6}, {7, 8}, {1, 2}}}, + {name: "multipoint + multipoint", g: mp1, h: mp2, want: geom.MultiPoint{{5, 6}, {7, 8}, {9, 9}}}, + + {name: "linestring + linestring", g: l1, h: l2, want: geom.MultiLineString{l1, l2}}, + { + name: "linestring + multilinestring", g: l1, h: ml1, + want: geom.MultiLineString{{{0, 0}, {1, 1}}, {{4, 4}, {5, 5}}}, + }, + { + name: "multilinestring + linestring", g: ml1, h: l1, + want: geom.MultiLineString{{{4, 4}, {5, 5}}, {{0, 0}, {1, 1}}}, + }, + { + name: "multilinestring + multilinestring", g: ml1, h: ml2, + want: geom.MultiLineString{{{4, 4}, {5, 5}}, {{6, 6}, {7, 7}}}, + }, + + {name: "polygon + polygon", g: poly1, h: poly2, want: geom.MultiPolygon{poly1, poly2}}, + { + name: "polygon + multipolygon", g: poly1, h: mpoly1, + want: geom.MultiPolygon{{{{0, 0}, {1, 0}, {1, 1}}}, {{{4, 4}, {5, 4}, {5, 5}}}}, + }, + { + name: "multipolygon + polygon", g: mpoly1, h: poly1, + want: geom.MultiPolygon{{{{4, 4}, {5, 4}, {5, 5}}}, {{{0, 0}, {1, 0}, {1, 1}}}}, + }, + { + name: "multipolygon + multipolygon", g: mpoly1, h: mpoly2, + want: geom.MultiPolygon{{{{4, 4}, {5, 4}, {5, 5}}}, {{{6, 6}, {7, 6}, {7, 7}}}}, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + got := MergeGeometries(tt.g, tt.h) + assert.Equal(t, tt.want, got) + }) + } +} + +func TestMergeGeometriesPanicsOnMismatchedTypes(t *testing.T) { + tests := []struct { + name string + g, h geom.Geometry + }{ + {name: "point with linestring", g: geom.Point{1, 2}, h: geom.LineString{{0, 0}, {1, 1}}}, + {name: "linestring with polygon", g: geom.LineString{{0, 0}, {1, 1}}, h: geom.Polygon{{{0, 0}, {1, 0}, {1, 1}}}}, + {name: "polygon with point", g: geom.Polygon{{{0, 0}, {1, 0}, {1, 1}}}, h: geom.Point{1, 2}}, + {name: "unknown geometry types", g: 42, h: 7}, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + assert.Panics(t, func() { MergeGeometries(tt.g, tt.h) }) + }) + } +} + +func TestGeometrySliceToGeom(t *testing.T) { + p1 := geom.Point{1, 2} + p2 := geom.Point{3, 4} + l1 := geom.LineString{{0, 0}, {1, 1}} + l2 := geom.LineString{{2, 2}, {3, 3}} + poly1 := geom.Polygon{{{0, 0}, {1, 0}, {1, 1}}} + poly2 := geom.Polygon{{{2, 2}, {3, 2}, {3, 3}}} + + tests := []struct { + name string + geometries []geom.Geometry + want geom.Geometry + }{ + { + name: "empty slice", + geometries: []geom.Geometry{}, + want: nil, + }, + { + name: "nil slice", + geometries: nil, + want: nil, + }, + { + name: "single point", + geometries: []geom.Geometry{p1}, + want: p1, + }, + { + name: "single linestring", + geometries: []geom.Geometry{l1}, + want: l1, + }, + { + name: "single polygon", + geometries: []geom.Geometry{poly1}, + want: poly1, + }, + { + name: "multiple points", + geometries: []geom.Geometry{p1, p2}, + want: geom.MultiPoint{p1, p2}, + }, + { + name: "multiple linestrings", + geometries: []geom.Geometry{l1, l2}, + want: geom.MultiLineString{l1, l2}, + }, + { + name: "multiple polygons", + geometries: []geom.Geometry{poly1, poly2}, + want: geom.MultiPolygon{poly1, poly2}, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + got := GeometrySliceToGeom(tt.geometries) + assert.Equal(t, tt.want, got) + }) + } +} + +func TestGeometrySliceToGeomPanics(t *testing.T) { + tests := []struct { + name string + geometries []geom.Geometry + }{ + { + name: "mixed points and linestrings", + geometries: []geom.Geometry{geom.Point{1, 2}, geom.LineString{{0, 0}, {1, 1}}}, + }, + { + name: "mixed linestrings and polygons", + geometries: []geom.Geometry{geom.LineString{{0, 0}, {1, 1}}, geom.Polygon{{{0, 0}, {1, 0}, {1, 1}}}}, + }, + { + name: "mixed points and polygons", + geometries: []geom.Geometry{geom.Point{1, 2}, geom.Polygon{{{0, 0}, {1, 0}, {1, 1}}}}, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + assert.Panics(t, func() { GeometrySliceToGeom(tt.geometries) }) + }) + } +} + +func TestPointSlice(t *testing.T) { + tests := []struct { + name string + g geom.Geometry + want []geom.Point + }{ + { + name: "point", + g: geom.Point{1, 2}, + want: []geom.Point{{1, 2}}, + }, + { + name: "line string", + g: geom.LineString{{0, 0}, {1, 1}, {2, 2}}, + want: []geom.Point{{0, 0}, {1, 1}, {2, 2}}, + }, + { + name: "polygon with single ring", + g: geom.Polygon{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + want: []geom.Point{{0, 0}, {1, 0}, {1, 1}, {0, 1}}, + }, + { + name: "polygon with hole", + g: geom.Polygon{ + {{0, 0}, {10, 0}, {10, 10}, {0, 10}}, + {{2, 2}, {4, 2}, {4, 4}, {2, 4}}, + }, + want: []geom.Point{ + {0, 0}, + {10, 0}, + {10, 10}, + {0, 10}, + {2, 2}, + {4, 2}, + {4, 4}, + {2, 4}, + }, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + got := PointSlice(tt.g) + assert.Equal(t, tt.want, got) + }) + } +} diff --git a/go.mod b/go.mod index 89a2d74..671277a 100644 --- a/go.mod +++ b/go.mod @@ -3,40 +3,56 @@ module github.com/pdok/texel go 1.26.2 require ( + github.com/Azure/azure-sdk-for-go/sdk/storage/azblob v1.8.1 + github.com/BurntSushi/toml v1.5.0 github.com/carlmjohnson/versioninfo v0.22.5 github.com/creasty/defaults v1.8.0 github.com/go-playground/validator/v10 v10.30.3 - github.com/go-spatial/geom v0.0.0-20220918193402-3cd2f5a9a082 + github.com/go-spatial/geom v0.1.0 + github.com/golang/protobuf v1.5.3 github.com/iancoleman/strcase v0.3.0 github.com/muesli/reflow v0.3.0 github.com/perimeterx/marshmallow v1.1.5 - github.com/stretchr/testify v1.11.1 + github.com/stretchr/testify v1.12.1 github.com/tobshub/go-sortedmap v1.0.3 github.com/urfave/cli/v2 v2.27.7 github.com/wk8/go-ordered-map/v2 v2.1.8 - golang.org/x/exp v0.0.0-20260611194520-c48552f49976 + golang.org/x/exp v0.0.0-20260813180055-c1d0aacb2297 ) require ( + github.com/Azure/azure-sdk-for-go/sdk/azcore v1.23.1 // indirect + github.com/Azure/azure-sdk-for-go/sdk/internal v1.12.0 // indirect + github.com/apache/arrow-go/v18 v18.7.0 // indirect + github.com/arolek/p v0.0.0-20191103215535-df3c295ed582 // indirect github.com/bahlo/generic-list-go v0.2.0 // indirect github.com/buger/jsonparser v1.1.1 // indirect + github.com/clipperhouse/uax29/v2 v2.7.0 // indirect github.com/cpuguy83/go-md2man/v2 v2.0.7 // indirect - github.com/davecgh/go-spew v1.1.1 // indirect github.com/gabriel-vasile/mimetype v1.4.13 // indirect github.com/gdey/errors v0.0.0-20190426172550-8ebd5bc891fb // indirect github.com/go-playground/locales v0.14.1 // indirect github.com/go-playground/universal-translator v0.18.1 // indirect + github.com/go-spatial/proj v0.3.0 // indirect + github.com/go-test/deep v1.1.1 // indirect + github.com/goccy/go-json v0.10.6 // indirect + github.com/google/flatbuffers v25.12.19+incompatible // indirect github.com/josharian/intern v1.0.0 // indirect + github.com/klauspost/compress v1.19.2 // indirect + github.com/klauspost/cpuid/v2 v2.4.0 // indirect github.com/leodido/go-urn v1.4.0 // indirect github.com/mailru/easyjson v0.7.7 // indirect - github.com/mattn/go-runewidth v0.0.12 // indirect - github.com/mattn/go-sqlite3 v1.14.17 // indirect - github.com/pmezard/go-difflib v1.0.0 // indirect - github.com/rivo/uniseg v0.2.0 // indirect + github.com/mattn/go-runewidth v0.0.20 // indirect + github.com/mattn/go-sqlite3 v1.14.22 // indirect + github.com/pierrec/lz4/v4 v4.1.28 // indirect github.com/russross/blackfriday/v2 v2.1.0 // indirect github.com/xrash/smetrics v0.0.0-20240521201337-686a1a2994c1 // indirect - golang.org/x/crypto v0.52.0 // indirect - golang.org/x/sys v0.45.0 // indirect - golang.org/x/text v0.37.0 // indirect + github.com/zeebo/xxh3 v1.1.0 // indirect + go.yaml.in/yaml/v3 v3.0.5 // indirect + golang.org/x/crypto v0.55.0 // indirect + golang.org/x/net v0.58.0 // indirect + golang.org/x/sys v0.47.0 // indirect + golang.org/x/text v0.41.0 // indirect + google.golang.org/protobuf v1.36.11 // indirect gopkg.in/yaml.v3 v3.0.1 // indirect ) diff --git a/go.sum b/go.sum index 4d8f73b..d07a0ed 100644 --- a/go.sum +++ b/go.sum @@ -1,3 +1,24 @@ +github.com/Azure/azure-sdk-for-go/sdk/azcore v1.23.1 h1:zvXfGJCWvywnCA814d8ZiVyt+fm9nnTE8xSb99zRyfo= +github.com/Azure/azure-sdk-for-go/sdk/azcore v1.23.1/go.mod h1:iptorS+VYKFL2N6PnebpS91dubG35eAOEERnT4PJbQU= +github.com/Azure/azure-sdk-for-go/sdk/azidentity v1.14.0 h1:CU4+EJeJi3TKYWEcYuSdWsjzw0nVsK/H0MSQOiPcymU= +github.com/Azure/azure-sdk-for-go/sdk/azidentity v1.14.0/go.mod h1:q0+UTSRvShwUCrR/s5HtyInYphN7Wvxb7snFM3u+SLA= +github.com/Azure/azure-sdk-for-go/sdk/internal v1.12.0 h1:fhqpLE3UEXi9lPaBRpQ6XuRW0nU7hgg4zlmZZa+a9q4= +github.com/Azure/azure-sdk-for-go/sdk/internal v1.12.0/go.mod h1:7dCRMLwisfRH3dBupKeNCioWYUZ4SS09Z14H+7i8ZoY= +github.com/Azure/azure-sdk-for-go/sdk/resourcemanager/storage/armstorage v1.8.1 h1:/Zt+cDPnpC3OVDm/JKLOs7M2DKmLRIIp3XIx9pHHiig= +github.com/Azure/azure-sdk-for-go/sdk/resourcemanager/storage/armstorage v1.8.1/go.mod h1:Ng3urmn6dYe8gnbCMoHHVl5APYz2txho3koEkV2o2HA= +github.com/Azure/azure-sdk-for-go/sdk/storage/azblob v1.8.1 h1:gkBLVmB3Z/HnGP/Jo4o12/RDpi0agnKav6sCKsX5Vu0= +github.com/Azure/azure-sdk-for-go/sdk/storage/azblob v1.8.1/go.mod h1:e3/1P5K+jIUi9JevDRklq/tFeTvbBb75bNAjU4xd31w= +github.com/AzureAD/microsoft-authentication-library-for-go v1.8.0 h1:Nljr4q1GRA/5vCrMONS+g4u4LRHNgOXVSh3O43J2CnI= +github.com/AzureAD/microsoft-authentication-library-for-go v1.8.0/go.mod h1:Y33QHnf0FfdVewFFISOGe20mkZbxX4H839o955/PoeI= +github.com/BurntSushi/toml v1.5.0 h1:W5quZX/G/csjUnuI8SUYlsHs9M38FC7znL0lIO+DvMg= +github.com/BurntSushi/toml v1.5.0/go.mod h1:ukJfTF/6rtPPRCnwkur4qwRxa8vTRFBF0uk2lLoLwho= +github.com/andybalholm/brotli v1.2.2 h1:HzTuoo2ErYQqf5qvcJInB8uvqSVxRttzkFexPWtnceM= +github.com/andybalholm/brotli v1.2.2/go.mod h1:rzTDkvFWvIrjDXZHkuS16NPggd91W3kUSvPlQ1pLaKY= +github.com/apache/arrow-go/v18 v18.7.0 h1:Vw/i+cJyebUofT7JlqFpe65LrmwxULn166jjwStM4HY= +github.com/apache/arrow-go/v18 v18.7.0/go.mod h1:PM6IigLJkdMwIpeHXnymo+xZ52f42a9EYiLtRel4p/A= +github.com/apache/thrift v0.24.0 h1:zy31L1a49QTNB2bG1BBfMXol3yJrTH975G3pPubQVLQ= +github.com/apache/thrift v0.24.0/go.mod h1:zPt6WxgvTOM6hF92y8C+MkEM5LMxZuk4JcQOiU4Esvs= +github.com/arolek/p v0.0.0-20191103215535-df3c295ed582 h1:DugKk4B3PpqfZ1QunDSPqNZDAErQmdCgwurloi1Axaw= github.com/arolek/p v0.0.0-20191103215535-df3c295ed582/go.mod h1:JPNItmi3yb44Q5QWM+Kh5n9oeRhfcJzPNS90mbLo25U= github.com/bahlo/generic-list-go v0.2.0 h1:5sz/EEAK+ls5wF+NeqDpk5+iNdMDXrh3z3nPnH1Wvgk= github.com/bahlo/generic-list-go v0.2.0/go.mod h1:2KvAjgMlE5NNynlg/5iLrrCCZ2+5xWbdbCW3pNTGyYg= @@ -5,13 +26,12 @@ github.com/buger/jsonparser v1.1.1 h1:2PnMjfWD7wBILjqQbt530v576A/cAbQvEW9gGIpYMU github.com/buger/jsonparser v1.1.1/go.mod h1:6RYKKt7H4d4+iWqouImQ9R2FZql3VbhNgx27UK13J/0= github.com/carlmjohnson/versioninfo v0.22.5 h1:O00sjOLUAFxYQjlN/bzYTuZiS0y6fWDQjMRvwtKgwwc= github.com/carlmjohnson/versioninfo v0.22.5/go.mod h1:QT9mph3wcVfISUKd0i9sZfVrPviHuSF+cUtLjm2WSf8= +github.com/clipperhouse/uax29/v2 v2.7.0 h1:+gs4oBZ2gPfVrKPthwbMzWZDaAFPGYK72F0NJv2v7Vk= +github.com/clipperhouse/uax29/v2 v2.7.0/go.mod h1:EFJ2TJMRUaplDxHKj1qAEhCtQPW2tJSwu5BF98AuoVM= github.com/cpuguy83/go-md2man/v2 v2.0.7 h1:zbFlGlXEAKlwXpmvle3d8Oe3YnkKIK4xSRTd3sHPnBo= github.com/cpuguy83/go-md2man/v2 v2.0.7/go.mod h1:oOW0eioCTA6cOiMLiUPZOpcVxMig6NIQQ7OS05n1F4g= github.com/creasty/defaults v1.8.0 h1:z27FJxCAa0JKt3utc0sCImAEb+spPucmKoOdLHvHYKk= github.com/creasty/defaults v1.8.0/go.mod h1:iGzKe6pbEHnpMPtfDXZEr0NVxWnPTjb1bbDy08fPzYM= -github.com/davecgh/go-spew v1.1.0/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38= -github.com/davecgh/go-spew v1.1.1 h1:vj9j/u1bqnvCEfJOwUhtlOARqs3+rkHYY13jYWTU97c= -github.com/davecgh/go-spew v1.1.1/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38= github.com/gabriel-vasile/mimetype v1.4.13 h1:46nXokslUBsAJE/wMsp5gtO500a4F3Nkz9Ufpk2AcUM= github.com/gabriel-vasile/mimetype v1.4.13/go.mod h1:d+9Oxyo1wTzWdyVUPMmXFvp4F9tea18J8ufA774AB3s= github.com/gdey/errors v0.0.0-20190426172550-8ebd5bc891fb h1:FYO+lZtAUnakgSW9xYs7QvgawjCDM5wgHaXoDhYHNH4= @@ -24,44 +44,65 @@ github.com/go-playground/universal-translator v0.18.1 h1:Bcnm0ZwsGyWbCzImXv+pAJn github.com/go-playground/universal-translator v0.18.1/go.mod h1:xekY+UJKNuX9WP91TpwSH2VMlDf28Uj24BCp08ZFTUY= github.com/go-playground/validator/v10 v10.30.3 h1:4MU6YkEwx7GbcPJOZxrtbu+QfF3pJLJuaYTeAH0DYy8= github.com/go-playground/validator/v10 v10.30.3/go.mod h1:4Axh7oCNGcoGkqLoE4YWt6n20mcEIsPRlB7vPk3lpyc= -github.com/go-spatial/geom v0.0.0-20220918193402-3cd2f5a9a082 h1:3+Swhq7I1stScm1DE4PUWKjGJRbi+VPvisj6tFz0prs= -github.com/go-spatial/geom v0.0.0-20220918193402-3cd2f5a9a082/go.mod h1:YU06tBGGstQCUX7vMuyF44RRneTdjGxOrp8OJHu4t9Q= -github.com/go-spatial/proj v0.2.0 h1:sii5Now3GFEyR9hV2SBJFWFz95s4xSaZmP0aYDk1K6c= -github.com/go-spatial/proj v0.2.0/go.mod h1:ePHHp7ITVc4eIVW5sgG/0Eu9RMAGOQUgM/D1ZkccY+0= -github.com/go-test/deep v1.0.8 h1:TDsG77qcSprGbC6vTN8OuXp5g+J+b5Pcguhf7Zt61VM= -github.com/go-test/deep v1.0.8/go.mod h1:5C2ZWiW0ErCdrYzpqxLbTX7MG14M9iiw8DgHncVwcsE= -github.com/golang/protobuf v1.3.2/go.mod h1:6lQm79b+lXiMfvg/cZm0SGofjICqVBUtrP5yJMmIC1U= -github.com/google/uuid v1.0.0/go.mod h1:TIyPZe4MgqvfeYDBFedMoGGpEw/LqOeaOT+nhxU+yHo= +github.com/go-spatial/geom v0.1.0 h1:IZ6LVz0Kkgd8+U83MmI0J4L4TBDzk5I8xiYe9Ci+aHk= +github.com/go-spatial/geom v0.1.0/go.mod h1:yNr22dnkQyFEwq2MJVmmiReiFesQpXKlosAG+qUSPus= +github.com/go-spatial/proj v0.3.0 h1:2IRb4elpnmVy1tUKRGBR68uiyZU4Pmny8yZC7TEUdSg= +github.com/go-spatial/proj v0.3.0/go.mod h1:wF0V0A60cTD5uzXeg+P52E1GtMY40PIiyYNJVvHXm6k= +github.com/go-test/deep v1.1.1 h1:0r/53hagsehfO4bzD2Pgr/+RgHqhmf+k1Bpse2cTu1U= +github.com/go-test/deep v1.1.1/go.mod h1:5C2ZWiW0ErCdrYzpqxLbTX7MG14M9iiw8DgHncVwcsE= +github.com/goccy/go-json v0.10.6 h1:p8HrPJzOakx/mn/bQtjgNjdTcN+/S6FcG2CTtQOrHVU= +github.com/goccy/go-json v0.10.6/go.mod h1:oq7eo15ShAhp70Anwd5lgX2pLfOS3QCiwU/PULtXL6M= +github.com/golang-jwt/jwt/v5 v5.3.1 h1:kYf81DTWFe7t+1VvL7eS+jKFVWaUnK9cB1qbwn63YCY= +github.com/golang-jwt/jwt/v5 v5.3.1/go.mod h1:fxCRLWMO43lRc8nhHWY6LGqRcf+1gQWArsqaEUEa5bE= +github.com/golang/protobuf v1.5.0/go.mod h1:FsONVRAS9T7sI+LIUmWTfcYkHO4aIWwzhcaSAoJOfIk= +github.com/golang/protobuf v1.5.3 h1:KhyjKVUg7Usr/dYsdSqoFveMYd5ko72D+zANwlG1mmg= +github.com/golang/protobuf v1.5.3/go.mod h1:XVQd3VNwM+JqD3oG2Ue2ip4fOMUkwXdXDdiuN0vRsmY= +github.com/google/flatbuffers v25.12.19+incompatible h1:haMV2JRRJCe1998HeW/p0X9UaMTK6SDo0ffLn2+DbLs= +github.com/google/flatbuffers v25.12.19+incompatible/go.mod h1:1AeVuKshWv4vARoZatz6mlQ0JxURH0Kv5+zNeJKJCa8= +github.com/google/go-cmp v0.5.5/go.mod h1:v8dTdLbMG2kIc/vJvl+f65V22dbkXbowE6jgT/gNBxE= +github.com/google/go-cmp v0.7.0 h1:wk8382ETsv4JYUZwIsn6YpYiWiBsYLSJiTsyBybVuN8= +github.com/google/go-cmp v0.7.0/go.mod h1:pXiqmnSA92OHEEa9HXL2W4E7lf9JzCmGVUdgjX3N/iU= +github.com/google/uuid v1.6.0 h1:NIvaJDMOsjHA8n1jAhLSgzrAzy1Hgr+hNrb57e+94F0= +github.com/google/uuid v1.6.0/go.mod h1:TIyPZe4MgqvfeYDBFedMoGGpEw/LqOeaOT+nhxU+yHo= github.com/iancoleman/strcase v0.3.0 h1:nTXanmYxhfFAMjZL34Ov6gkzEsSJZ5DbhxWjvSASxEI= github.com/iancoleman/strcase v0.3.0/go.mod h1:iwCmte+B7n89clKwxIoIXy/HfoL7AsD47ZCWhYzw7ho= github.com/josharian/intern v1.0.0 h1:vlS4z54oSdjm0bgjRigI+G1HpF+tI+9rE5LLzOg8HmY= github.com/josharian/intern v1.0.0/go.mod h1:5DoeVV0s6jJacbCEi61lwdGj/aVlrQvzHFFd8Hwg//Y= +github.com/klauspost/compress v1.19.2 h1:hMRETovs/pu/dVWN7zIT1PGG8t509MwT6bO7XSi26R8= +github.com/klauspost/compress v1.19.2/go.mod h1:cwPg85FWrGar70rWktvGQj8/hthj3wpl0PGDogxkrSQ= +github.com/klauspost/cpuid/v2 v2.4.0 h1:S6Hrbc7+ywsr0r+RLapfGBHfyefhCTwEh3A0tV913Dw= +github.com/klauspost/cpuid/v2 v2.4.0/go.mod h1:19jmZ9mjzoF//ddRSUsv0zfBTJWh3QJh9FNxZTMrGxU= +github.com/kr/pretty v0.3.1 h1:flRD4NNwYAUpkphVc1HcthR4KEIFJ65n8Mw5qdRn3LE= +github.com/kr/pretty v0.3.1/go.mod h1:hoEshYVHaxMs3cyo3Yncou5ZscifuDolrwPKZanG3xk= +github.com/kr/text v0.2.0 h1:5Nx0Ya0ZqY2ygV366QzturHI13Jq95ApcVaJBhpS+AY= +github.com/kr/text v0.2.0/go.mod h1:eLer722TekiGuMkidMxC/pM04lWEeraHUUmBw8l2grE= +github.com/kylelemons/godebug v1.1.0 h1:RPNrshWIDI6G2gRW9EHilWtl7Z6Sb1BR0xunSBf0SNc= +github.com/kylelemons/godebug v1.1.0/go.mod h1:9/0rRGxNHcop5bhtWyNeEfOS8JIWk580+fNqagV/RAw= github.com/leodido/go-urn v1.4.0 h1:WT9HwE9SGECu3lg4d/dIA+jxlljEa1/ffXKmRjqdmIQ= github.com/leodido/go-urn v1.4.0/go.mod h1:bvxc+MVxLKB4z00jd1z+Dvzr47oO32F/QSNjSBOlFxI= github.com/mailru/easyjson v0.7.7 h1:UGYAvKxe3sBsEDzO8ZeWOSlIQfWFlxbzLZe7hwFURr0= github.com/mailru/easyjson v0.7.7/go.mod h1:xzfreul335JAWq5oZzymOObrkdz5UnU4kGfJJLY9Nlc= -github.com/mattn/go-runewidth v0.0.12 h1:Y41i/hVW3Pgwr8gV+J23B9YEY0zxjptBuCWEaxmAOow= github.com/mattn/go-runewidth v0.0.12/go.mod h1:RAqKPSqVFrSLVXbA8x7dzmKdmGzieGRCM46jaSJTDAk= -github.com/mattn/go-sqlite3 v1.12.0/go.mod h1:FPy6KqzDD04eiIsT53CuJW3U88zkxoIYsOqkbpncsNc= -github.com/mattn/go-sqlite3 v1.14.17 h1:mCRHCLDUBXgpKAqIKsaAaAsrAlbkeomtRFKXh2L6YIM= -github.com/mattn/go-sqlite3 v1.14.17/go.mod h1:2eHXhiwb8IkHr+BDWZGa96P6+rkvnG63S2DGjv9HUNg= -github.com/mattn/goveralls v0.0.3-0.20180319021929-1c14a4061c1c/go.mod h1:8d1ZMHsd7fW6IRPKQh46F2WRpyib5/X4FOpevwGNQEw= +github.com/mattn/go-runewidth v0.0.20 h1:WcT52H91ZUAwy8+HUkdM3THM6gXqXuLJi9O3rjcQQaQ= +github.com/mattn/go-runewidth v0.0.20/go.mod h1:XBkDxAl56ILZc9knddidhrOlY5R/pDhgLpndooCuJAs= +github.com/mattn/go-sqlite3 v1.14.22 h1:2gZY6PC6kBnID23Tichd1K+Z0oS6nE/XwU+Vz/5o4kU= +github.com/mattn/go-sqlite3 v1.14.22/go.mod h1:Uh1q+B4BYcTPb+yiD3kU8Ct7aC0hY9fxUwlHK0RXw+Y= github.com/muesli/reflow v0.3.0 h1:IFsN6K9NfGtjeggFP+68I4chLZV2yIKsXJFNZ+eWh6s= github.com/muesli/reflow v0.3.0/go.mod h1:pbwTDkVPibjO2kyvBQRBxTWEEGDGq0FlB1BIKtnHY/8= -github.com/pborman/uuid v1.2.0/go.mod h1:X/NO0urCmaxf9VXbdlT7C2Yzkj2IKimNn4k+gtPdI/k= github.com/perimeterx/marshmallow v1.1.5 h1:a2LALqQ1BlHM8PZblsDdidgv1mWi1DgC2UmX50IvK2s= github.com/perimeterx/marshmallow v1.1.5/go.mod h1:dsXbUu8CRzfYP5a87xpp0xq9S3u0Vchtcl8we9tYaXw= -github.com/pmezard/go-difflib v1.0.0 h1:4DBwDE0NGyQoBHbLQYPwSUPoCMWR5BEzIk/f1lZbAQM= -github.com/pmezard/go-difflib v1.0.0/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZNVY4sRDYZ/4= +github.com/pierrec/lz4/v4 v4.1.28 h1:pPEPwRJ4kybBTfGt28q7lQsRJQHhC08axprdLD5Ppio= +github.com/pierrec/lz4/v4 v4.1.28/go.mod h1:EoQMVJgeeEOMsCqCzqFm2O0cJvljX2nGZjcRIPL34O4= +github.com/pkg/browser v0.0.0-20240102092130-5ac0b6a4141c h1:+mdjkGKdHQG3305AYmdv1U2eRNDiU2ErMBj1gwrq8eQ= +github.com/pkg/browser v0.0.0-20240102092130-5ac0b6a4141c/go.mod h1:7rwL4CYBLnjLxUqIJNnCWiEdr3bn6IUYi15bNlnbCCU= github.com/rivo/uniseg v0.1.0/go.mod h1:J6wj4VEh+S6ZtnVlnTBMWIodfgj8LQOQFoIToxlJtxc= -github.com/rivo/uniseg v0.2.0 h1:S1pD9weZBuJdFmowNwbpi7BJ8TNftyUImj/0WQi72jY= github.com/rivo/uniseg v0.2.0/go.mod h1:J6wj4VEh+S6ZtnVlnTBMWIodfgj8LQOQFoIToxlJtxc= +github.com/rogpeppe/go-internal v1.12.0 h1:exVL4IDcn6na9z1rAb56Vxr+CgyK3nn3O+epU5NdKM8= +github.com/rogpeppe/go-internal v1.12.0/go.mod h1:E+RYuTGaKKdloAfM02xzb0FW3Paa99yedzYV+kq4uf4= github.com/russross/blackfriday/v2 v2.1.0 h1:JIOH55/0cWyOuilr9/qlrm0BSXldqnqwMsf35Ld67mk= github.com/russross/blackfriday/v2 v2.1.0/go.mod h1:+Rmxgy9KzJVeS9/2gXHxylqXiyQDYRxCVz55jmeOWTM= -github.com/stretchr/objx v0.1.0/go.mod h1:HFkY916IF+rwdDfMAkV7OtwuqBVzrE8GR6GFx+wExME= -github.com/stretchr/testify v1.4.0/go.mod h1:j7eGeouHqKxXV5pUuKE4zz7dFj8WfuZ+81PSLYec5m4= -github.com/stretchr/testify v1.11.1 h1:7s2iGBzp5EwR7/aIZr8ao5+dra3wiQyKjjFuvgVKu7U= -github.com/stretchr/testify v1.11.1/go.mod h1:wZwfW3scLgRK+23gO65QZefKpKQRnfz6sD981Nm4B6U= +github.com/stretchr/testify v1.12.1 h1:EuwCh5fleGS7H32xRwO3wRGT7DxrDhLAT6FF8MpWDWE= +github.com/stretchr/testify v1.12.1/go.mod h1:MDEgiDPPsNp5cuIrHPPCyornHKgEVbtFUmoNlxoYthg= github.com/tobshub/go-sortedmap v1.0.3 h1:oUhj/5tqzjTX4bhWqB1ZFTDtMULJ1ZYUnS8WAugSfjY= github.com/tobshub/go-sortedmap v1.0.3/go.mod h1:JLxyU94+lfKuCgelxXpwRr29ei6SqLbaeVuNVMvENbE= github.com/ugorji/go/codec v1.2.7 h1:YPXUKf7fYbp/y8xloBqZOw2qaVggbfwMlI8WM3wZUJ0= @@ -72,23 +113,33 @@ github.com/wk8/go-ordered-map/v2 v2.1.8 h1:5h/BUHu93oj4gIdvHHHGsScSTMijfx5PeYkE/ github.com/wk8/go-ordered-map/v2 v2.1.8/go.mod h1:5nJHM5DyteebpVlHnWMV0rPz6Zp7+xBAnxjb1X5vnTw= github.com/xrash/smetrics v0.0.0-20240521201337-686a1a2994c1 h1:gEOO8jv9F4OT7lGCjxCBTO/36wtF6j2nSip77qHd4x4= github.com/xrash/smetrics v0.0.0-20240521201337-686a1a2994c1/go.mod h1:Ohn+xnUBiLI6FVj/9LpzZWtj1/D6lUovWYBkxHVV3aM= -golang.org/x/crypto v0.0.0-20190308221718-c2843e01d9a2/go.mod h1:djNgcEr1/C05ACkg1iLfiJU5Ep61QUkGW8qpdssI0+w= -golang.org/x/crypto v0.52.0 h1:RMs7fP2rXdep0CftQlK8Uf+kibLm7qkCcradZWYz988= -golang.org/x/crypto v0.52.0/go.mod h1:1QgfPxDqh0T2M/elOJtp9RvuR95kVjir0e6/BvEmGbc= -golang.org/x/exp v0.0.0-20260611194520-c48552f49976 h1:X8Hz2ImujgbmetVuW+w2YkyZChE3cBpZi2P158rTG9M= -golang.org/x/exp v0.0.0-20260611194520-c48552f49976/go.mod h1:vnf4pv9iKZXY58sQE1L86zmNWJ4159e1RkcWiLCkeEY= -golang.org/x/net v0.0.0-20190620200207-3b0461eec859/go.mod h1:z5CRVTTTmAJ677TzLLGU+0bjPO0LkuOLi4/5GtJWs/s= -golang.org/x/sync v0.0.0-20190423024810-112230192c58/go.mod h1:RxMgew5VJxzue5/jJTE5uejpjVlOe/izrB70Jof72aM= -golang.org/x/sys v0.0.0-20190215142949-d0b11bdaac8a/go.mod h1:STP8DvDyc/dI5b8T5hshtkjS+E42TnysNCUPdjciGhY= -golang.org/x/sys v0.45.0 h1:dO4czNzziLiiXplLQgBCEpCvXQ3dnkn0SdaZSYdQ+FY= -golang.org/x/sys v0.45.0/go.mod h1:4GL1E5IUh+htKOUEOaiffhrAeqysfVGipDYzABqnCmw= -golang.org/x/text v0.3.0/go.mod h1:NqM8EUOU14njkJ3fqMW+pc6Ldnwhi/IjpwHt7yyuwOQ= -golang.org/x/text v0.37.0 h1:Cqjiwd9eSg8e0QAkyCaQTNHFIIzWtidPahFWR83rTrc= -golang.org/x/text v0.37.0/go.mod h1:a5sjxXGs9hsn/AJVwuElvCAo9v8QYLzvavO5z2PiM38= -golang.org/x/tools v0.0.0-20191114222411-4191b8cbba09/go.mod h1:b+2E5dAYhXwXZwtnZ6UAqBI28+e2cm9otk0dWdXHAEo= -golang.org/x/xerrors v0.0.0-20190717185122-a985d3407aa7/go.mod h1:I/5z698sn9Ka8TeJc9MKroUUfqBBauWjQqLJ2OPfmY0= -gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405 h1:yhCVgyC4o1eVCa2tZl7eS0r+SDo693bJlVdllGtEeKM= +github.com/zeebo/assert v1.3.0 h1:g7C04CbJuIDKNPFHmsk4hwZDO5O+kntRxzaUoNXj+IQ= +github.com/zeebo/assert v1.3.0/go.mod h1:Pq9JiuJQpG8JLJdtkwrJESF0Foym2/D9XMU5ciN/wJ0= +github.com/zeebo/xxh3 v1.1.0 h1:s7DLGDK45Dyfg7++yxI0khrfwq9661w9EN78eP/UZVs= +github.com/zeebo/xxh3 v1.1.0/go.mod h1:IisAie1LELR4xhVinxWS5+zf1lA4p0MW4T+w+W07F5s= +go.yaml.in/yaml/v3 v3.0.5 h1:N6y/pJk8buWs9NY5ERU2HSMfm+IuD/OtfdAnq6kESPw= +go.yaml.in/yaml/v3 v3.0.5/go.mod h1:HVTZu1O7/Vkt2N+BFy8Zza+lnLsABggaTM2ZpNIGuKg= +golang.org/x/crypto v0.55.0 h1:+KWHjbgOaAQ66dh/YlkZKHlz9ZUlq61AFirAR9ntP8M= +golang.org/x/crypto v0.55.0/go.mod h1:uq0V9dE/fzQuJtbnL+2EhWOE63vo164FY8xqEnV9xis= +golang.org/x/exp v0.0.0-20260813180055-c1d0aacb2297 h1:YXnL44eJ77R+ji4/ooy8UsXIhz+lbi2Qgdlc8iRN0gY= +golang.org/x/exp v0.0.0-20260813180055-c1d0aacb2297/go.mod h1:Mkmymgv+uMpSQ/XxJ/7GpdrdYoqm3u72jEbpCLiJmNk= +golang.org/x/net v0.58.0 h1:ynWG7rqYi4ccpTEuPZ2QGWHktVEM9DMCj9yzDE0Q7To= +golang.org/x/net v0.58.0/go.mod h1:YwCddHnFlT7eLQqVprV19OnhLGtc5xOKgE0RyqgfWAU= +golang.org/x/sync v0.22.0 h1:SZjpbeLmrCk4xhRSZFNZW5gFUeCeFgjekvI/+gfScek= +golang.org/x/sync v0.22.0/go.mod h1:9xrNwdLfx4jkKbNva9FpL6vEN7evnE43NNNJQ2LF3+0= +golang.org/x/sys v0.47.0 h1:o7XGOvZQCADBQQ4Y7VNq2dRWQR7JmOUW8Kxx4ZsNgWs= +golang.org/x/sys v0.47.0/go.mod h1:4GL1E5IUh+htKOUEOaiffhrAeqysfVGipDYzABqnCmw= +golang.org/x/text v0.41.0 h1:vz/seA0lnX87Othu2f/0L24RcgrXD9/YFTSuGjj3rH8= +golang.org/x/text v0.41.0/go.mod h1:jvf1O8ajNzZqhSrQBPbutR/EB83Cc0CFrezNQIwbb5M= +golang.org/x/xerrors v0.0.0-20191204190536-9bdfabe68543/go.mod h1:I/5z698sn9Ka8TeJc9MKroUUfqBBauWjQqLJ2OPfmY0= +gonum.org/v1/gonum v0.17.0 h1:VbpOemQlsSMrYmn7T2OUvQ4dqxQXU+ouZFQsZOx50z4= +gonum.org/v1/gonum v0.17.0/go.mod h1:El3tOrEuMpv2UdMrbNlKEh9vd86bmQ6vqIcDwxEOc1E= +google.golang.org/protobuf v1.26.0-rc.1/go.mod h1:jlhhOSvTdKEhbULTjvd4ARK9grFBp09yW+WbY/TyQbw= +google.golang.org/protobuf v1.26.0/go.mod h1:9q0QmTI4eRPtz6boOQmLYwt+qCgq0jsYwAQnmE0givc= +google.golang.org/protobuf v1.36.11 h1:fV6ZwhNocDyBLK0dj+fg8ektcVegBBuEolpbTQyBNVE= +google.golang.org/protobuf v1.36.11/go.mod h1:HTf+CrKn2C3g5S8VImy6tdcUvCska2kB7j23XfzDpco= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0= -gopkg.in/yaml.v2 v2.2.2/go.mod h1:hI93XBmqTisBFMUTm0b8Fm+jr3Dg1NNxqwp+5A1VGuI= +gopkg.in/check.v1 v1.0.0-20201130134442-10cb98267c6c h1:Hei/4ADfdWqJk1ZMxUNpqntNwaWcugrBjAiHlqqRiVk= +gopkg.in/check.v1 v1.0.0-20201130134442-10cb98267c6c/go.mod h1:JHkPIbrfpd72SG/EVd6muEfDQjcINNoR0C8j2r3qZ4Q= gopkg.in/yaml.v3 v3.0.1 h1:fxVm/GzAzEWqLHuvctI91KS9hhNmmWOoWu0XTYJS7CA= gopkg.in/yaml.v3 v3.0.1/go.mod h1:K4uyk7z7BCEPqu6E+C64Yfv1cQ7kz7rIZviUmN+EgEM= diff --git a/intgeom/intgeom.go b/intgeom/intgeom.go index 8285e28..00ac35e 100644 --- a/intgeom/intgeom.go +++ b/intgeom/intgeom.go @@ -56,7 +56,7 @@ func FromGeomOrd(o float64) M { // TODO implement this with integers func SegmentIntersect(l1, l2 Line) ([2]int64, bool) { intersection, intersects := planar.SegmentIntersect(l1.ToGeomLine(), l2.ToGeomLine()) - intIntersection := [2]int64{FromGeomOrd(intersection[0]), FromGeomOrd(intersection[0])} + intIntersection := [2]int64{FromGeomOrd(intersection[0]), FromGeomOrd(intersection[1])} return intIntersection, intersects } diff --git a/main.go b/main.go index 6b5c732..df434ec 100644 --- a/main.go +++ b/main.go @@ -10,10 +10,10 @@ import ( "slices" "syscall" + "github.com/pdok/texel/config" + "github.com/pdok/texel/mapslicehelp" "github.com/pdok/texel/pointindex" - "github.com/go-spatial/geom" - "github.com/carlmjohnson/versioninfo" "github.com/pdok/texel/processing" @@ -21,19 +21,29 @@ import ( "github.com/iancoleman/strcase" "github.com/pdok/texel/processing/gpkg" + "github.com/pdok/texel/processing/writetile" "github.com/pdok/texel/snap" "github.com/urfave/cli/v2" ) -const SOURCE string = `sourceGpkg` -const TARGET string = `targetGpkg` -const OVERWRITE string = `overwrite` -const TILEMATRIXSET string = `tilematrixset` -const TILEMATRICES string = `tilematrices` -const PAGESIZE string = `pagesize` -const KEEPPOINTSANDLINES string = `keeppointsandlines` -const IGNOREOUTSIDEGRID string = `ignoreoutsidegrid` -const REVERSEWINDINGORDER string = `reversewindingorder` +const ( + SOURCE string = `sourceGpkg` + TARGET string = `targetGpkg` + OVERWRITE string = `overwrite` + TILEMATRIXSET string = `tilematrixset` + TILEMATRICES string = `tilematrices` + PAGESIZE string = `pagesize` + KEEPPOINTSANDLINES string = `keeppointsandlines` + IGNOREOUTSIDEGRID string = `ignoreoutsidegrid` + REVERSEWINDINGORDER string = `reversewindingorder` + ENCODETILES string = `encodetiles` + TILEBUFFER string = `tilebuffer` + USELINETRACE string = `uselinetrace` + CLIP string = `clip` + + MVTCONFIG string = `config` + TILEMATRIX string = `tilematrix` +) //nolint:funlen func main() { @@ -42,148 +52,229 @@ func main() { app.Usage = "A Golang Polygon Snapping application" app.Version = versioninfo.Short() - app.Flags = []cli.Flag{ - &cli.StringFlag{ - Name: SOURCE, - Aliases: []string{"s"}, - Usage: "Source GPKG", - Required: true, - EnvVars: []string{strcase.ToScreamingSnake(SOURCE)}, - }, - &cli.StringFlag{ - Name: TARGET, - Aliases: []string{"t"}, - Usage: "Target GPKG (prefix). One GPKG per tile matrix cq zoom level will be created and the filename will be suffixed. E.g. target_6.gpkg", - Required: true, - EnvVars: []string{strcase.ToScreamingSnake(TARGET)}, - }, - &cli.BoolFlag{ - Name: OVERWRITE, - Aliases: []string{"o"}, - Usage: "Overwrite a target GPKG if it exists", - Required: false, - EnvVars: []string{strcase.ToScreamingSnake(OVERWRITE)}, - }, - &cli.StringFlag{ - Name: TILEMATRIXSET, - Aliases: []string{"tms"}, - Usage: `ID of a (built-in) tile matrix set. E.g.: NetherlandsRDNewQuad`, - Required: true, - EnvVars: []string{strcase.ToScreamingSnake(TILEMATRIXSET)}, - }, - &cli.StringFlag{ - Name: TILEMATRICES, - Aliases: []string{"z"}, - Usage: `IDs (usually the same as the zoom levels) of the tile matrices in the tile matrix set that should be processed for. JSON array of integers. E.g.: [4,5,6,7,8]`, - Required: true, - EnvVars: []string{strcase.ToScreamingSnake(TILEMATRICES)}, - }, - &cli.IntFlag{ - Name: PAGESIZE, - Aliases: []string{"p"}, - Usage: "Page Size, how many features are written per transaction to a target GPKG", - Value: 1000, - Required: false, - EnvVars: []string{strcase.ToScreamingSnake(PAGESIZE)}, - }, - &cli.BoolFlag{ - Name: KEEPPOINTSANDLINES, - Aliases: []string{"pl"}, - Usage: "Parts of polygons are reduced to points and lines after texel, keep these details or not.", - // TODO: This results in broken tiles, we should change this to a collapse option (optionally collapse rings from polygons to lines and points; grow lines as long as possible and add an option for line-length-to-keep) - Value: false, - Required: false, - EnvVars: []string{strcase.ToScreamingSnake(KEEPPOINTSANDLINES)}, - }, - &cli.BoolFlag{ - Name: IGNOREOUTSIDEGRID, - Aliases: []string{"iog"}, - Usage: "Ignore polygons that fall (partly) outside the grid, instead of panicking", - Value: false, - Required: false, - EnvVars: []string{strcase.ToScreamingSnake(IGNOREOUTSIDEGRID)}, + // `texel` without command defaults to `snap`. + app.DefaultCommand = "snap" + + // Define two commands: `texel snap` and `texel mvt`. + app.Commands = []*cli.Command{ + { + Name: "snap", + Usage: "Snap polygons in a source GeoPackage to a tile grid and write the result to a target GeoPackage", + Flags: []cli.Flag{ + &cli.StringFlag{ + Name: SOURCE, + Aliases: []string{"s"}, + Usage: "Source GPKG", + Required: true, + EnvVars: []string{strcase.ToScreamingSnake(SOURCE)}, + }, + &cli.StringFlag{ + Name: TARGET, + Aliases: []string{"t"}, + Usage: "Target GPKG (prefix). One GPKG per tile matrix cq zoom level will be created and the filename will be suffixed. E.g. target_6.gpkg", + Required: true, + EnvVars: []string{strcase.ToScreamingSnake(TARGET)}, + }, + &cli.BoolFlag{ + Name: OVERWRITE, + Aliases: []string{"o"}, + Usage: "Overwrite a target GPKG if it exists", + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(OVERWRITE)}, + }, + &cli.StringFlag{ + Name: TILEMATRIXSET, + Aliases: []string{"tms"}, + Usage: `ID of a (built-in) tile matrix set. E.g.: NetherlandsRDNewQuad`, + Required: true, + EnvVars: []string{strcase.ToScreamingSnake(TILEMATRIXSET)}, + }, + &cli.StringFlag{ + Name: TILEMATRICES, + Aliases: []string{"z"}, + Usage: `IDs (usually the same as the zoom levels) of the tile matrices in the tile matrix set that should be processed for. JSON array of integers. E.g.: [4,5,6,7,8]`, + Required: true, + EnvVars: []string{strcase.ToScreamingSnake(TILEMATRICES)}, + }, + &cli.IntFlag{ + Name: PAGESIZE, + Aliases: []string{"p"}, + Usage: "Page Size, how many features are written per transaction to a target GPKG", + Value: 1000, + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(PAGESIZE)}, + }, + &cli.BoolFlag{ + Name: KEEPPOINTSANDLINES, + Aliases: []string{"pl"}, + Usage: "Parts of polygons are reduced to points and lines after texel, keep these details or not.", + // TODO: This results in broken tiles, we should change this to a collapse option (optionally collapse rings from polygons to lines and points; grow lines as long as possible and add an option for line-length-to-keep) + Value: false, + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(KEEPPOINTSANDLINES)}, + }, + &cli.BoolFlag{ + Name: IGNOREOUTSIDEGRID, + Aliases: []string{"iog"}, + Usage: "Ignore polygons that fall (partly) outside the grid, instead of panicking", + Value: false, + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(IGNOREOUTSIDEGRID)}, + }, + &cli.BoolFlag{ + Name: REVERSEWINDINGORDER, + Aliases: []string{"rwo"}, + Usage: "Reverse the winding order of rings in polygons, contrary to the OGC simple features spec", + Value: false, + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(REVERSEWINDINGORDER)}, + }, + &cli.BoolFlag{ + Name: ENCODETILES, + Aliases: []string{"enc"}, + Usage: "Add tables with mapbox-encoded geometries.", + Value: false, + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(ENCODETILES)}, + }, + &cli.UintFlag{ + Name: TILEBUFFER, + Aliases: []string{"buf"}, + Usage: "Buffer (in internal pixels) used to select the tiles a polygon's geometry touches.", + Value: 0, + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(TILEBUFFER)}, + }, + &cli.BoolFlag{ + Name: USELINETRACE, + Aliases: []string{"lt"}, + Usage: "Use precise line-trace tile classification instead of a simple buffered bounding box to select tiles.", + Value: false, + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(USELINETRACE)}, + }, + &cli.StringFlag{ + Name: CLIP, + Aliases: []string{"c"}, + Usage: "Specify clip as \"[z,x,y]\", which denotes a quadrant in the target grid. Only encode geometries for tiles within this clip. Default: entire grid.", + Value: "[0,0,0]", + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(CLIP)}, + }, + }, + Action: runSnap, }, - &cli.BoolFlag{ - Name: REVERSEWINDINGORDER, - Aliases: []string{"rwo"}, - Usage: "Reverse the winding order of rings in polygons, contrary to the OGC simple features spec", - Value: false, - Required: false, - EnvVars: []string{strcase.ToScreamingSnake(REVERSEWINDINGORDER)}, + { + Name: "mvt", + Usage: "Build MVT tiles from an already-encoded target GeoPackage (see --" + ENCODETILES + ")", + Flags: []cli.Flag{ + &cli.StringFlag{ + Name: MVTCONFIG, + Aliases: []string{"c"}, + Usage: "Config toml file describing datasources, tilesets and layers", + Required: true, + EnvVars: []string{strcase.ToScreamingSnake(MVTCONFIG)}, + }, + &cli.UintFlag{ + Name: TILEMATRIX, + Aliases: []string{"z"}, + Usage: "Target zoomlevel", + Required: true, + EnvVars: []string{strcase.ToScreamingSnake(TILEMATRIX)}, + }, + }, + Action: func(c *cli.Context) error { + return runMVT(c.String(MVTCONFIG), c.Uint(TILEMATRIX)) + }, }, } - app.Action = func(c *cli.Context) error { - tileMatrixSet, err := tms20.LoadEmbeddedTileMatrixSet(c.String(TILEMATRIXSET)) - if err != nil { - return err - } - var tileMatrixIDs []int - err = json.Unmarshal([]byte(c.String(TILEMATRICES)), &tileMatrixIDs) - if err != nil { - return err - } - if err = validateTileMatrixSet(tileMatrixSet, tileMatrixIDs); err != nil { - return err - } + err := app.Run(os.Args) + if err != nil { + log.Fatal(err) + } +} - _, err = os.Stat(c.String(SOURCE)) - if os.IsNotExist(err) { - log.Fatalf("error opening source GeoPackage: %s", err) - } +func runSnap(c *cli.Context) error { + tileMatrixSet, err := tms20.LoadEmbeddedTileMatrixSet(c.String(TILEMATRIXSET)) + if err != nil { + return err + } + var tileMatrixIDs []int + err = json.Unmarshal([]byte(c.String(TILEMATRICES)), &tileMatrixIDs) + if err != nil { + return err + } + if err = validateTileMatrixSet(tileMatrixSet, tileMatrixIDs); err != nil { + return err + } - source := gpkg.SourceGeopackage{} - source.Init(c.String(SOURCE)) - defer source.Close() + _, err = os.Stat(c.String(SOURCE)) + if os.IsNotExist(err) { + log.Fatalf("error opening source GeoPackage: %s", err) + } + + clip, err := readClip(c.String(CLIP)) + if err != nil { + return err + } - targetPathFmt := injectSuffixIntoPath(c.String(TARGET)) + source := gpkg.SourceGeopackage{} + source.Init(c.String(SOURCE)) + defer source.Close() - gpkgTargets := make(map[int]*gpkg.TargetGeopackage, len(tileMatrixIDs)) - overwrite := c.Bool(OVERWRITE) - pagesize := c.Int(PAGESIZE) // TODO divide by tile matrices count - snapConfig := snap.Config{ - KeepPointsAndLines: c.Bool(KEEPPOINTSANDLINES), - IgnoreOutsideGrid: c.Bool(IGNOREOUTSIDEGRID), - ReverseWindingOrder: c.Bool(REVERSEWINDINGORDER), - } - for _, tmID := range tileMatrixIDs { - gpkgTargets[tmID] = initGPKGTarget(targetPathFmt, tmID, overwrite, pagesize) - defer gpkgTargets[tmID].Close() // yes, supposed to go here, want to close all at end of func - } + targetPathFmt := injectSuffixIntoPath(c.String(TARGET)) - tables := source.GetTableInfo() - for _, target := range gpkgTargets { - err = target.CreateTables(tables) - if err != nil { - log.Fatalf("error initialization the target GeoPackage: %s", err) - } - } + gpkgTargets := make(map[int]*gpkg.TargetGeopackage, len(tileMatrixIDs)) + overwrite := c.Bool(OVERWRITE) + pagesize := c.Int(PAGESIZE) // TODO divide by tile matrices count + snapConfig := processing.Config{ + KeepPointsAndLines: c.Bool(KEEPPOINTSANDLINES), + IgnoreOutsideGrid: c.Bool(IGNOREOUTSIDEGRID), + ReverseWindingOrder: c.Bool(REVERSEWINDINGORDER), + EncodeTiles: c.Bool(ENCODETILES), + Buffer: c.Uint(TILEBUFFER), + UseLineTrace: c.Bool(USELINETRACE), + Clip: processing.Clip{Z: clip[0], X: uint(clip[1]), Y: uint(clip[2])}, //nolint: gosec // G115 + } - log.Println("=== start snapping ===") + for _, tmID := range tileMatrixIDs { + gpkgTargets[tmID] = initGPKGTarget(targetPathFmt, tmID, overwrite, pagesize, c.Bool(ENCODETILES)) + defer gpkgTargets[tmID].Close() // yes, supposed to go here, want to close all at end of func + } - // need a copied map because of type difference processing.Target vs gpkg.TargetGeopackage - targets := make(map[int]processing.Target, len(gpkgTargets)) - for tmID, target := range gpkgTargets { - targets[tmID] = target - } - // Process the tables sequentially - for _, table := range tables { - log.Printf(" snapping %s", table.Name) - for _, target := range gpkgTargets { - source.Table = table - target.Table = table - } - processBySnapping(source, targets, tileMatrixSet, snapConfig) - log.Printf(" finished %s", table.Name) + tables := source.GetTableInfo() + for _, target := range gpkgTargets { + err = target.CreateTables(tables) + if err != nil { + return fmt.Errorf("error initialization the target GeoPackage: %w", err) } + } + + log.Println("=== start snapping ===") + processSnapTables(tables, source, gpkgTargets, snapConfig, tileMatrixSet) + log.Println("=== done snapping ===") - log.Println("=== done snapping ===") - return nil + return nil +} + +// Loop over tables, then loop over targets, and process +func processSnapTables(tables []gpkg.Table, source gpkg.SourceGeopackage, gpkgTargets map[int]*gpkg.TargetGeopackage, snapConfig processing.Config, tileMatrixSet tms20.TileMatrixSet) { + // need a copied map because of type difference processing.Target vs gpkg.TargetGeopackage + targets := make(map[int]processing.Target, len(gpkgTargets)) + for tmID, target := range gpkgTargets { + targets[tmID] = target } - err := app.Run(os.Args) - if err != nil { - log.Fatal(err) + // Process the tables sequentially + for _, table := range tables { + log.Printf(" snapping %s", table.Name) + for _, target := range gpkgTargets { + source.Table = table + target.Table = table + } + processBySnapping(source, targets, tileMatrixSet, snapConfig) + log.Printf(" finished %s", table.Name) } } @@ -200,7 +291,7 @@ func validateTileMatrixSet(tms tms20.TileMatrixSet, tileMatrixIDs []tms20.TMID) return pointindex.IsQuadTree(tms) } -func initGPKGTarget(targetPathFmt string, tmID int, overwrite bool, pagesize int) *gpkg.TargetGeopackage { +func initGPKGTarget(targetPathFmt string, tmID int, overwrite bool, pagesize int, encodeTiles bool) *gpkg.TargetGeopackage { targetPath := fmt.Sprintf(targetPathFmt, tmID) if overwrite { err := os.Remove(targetPath) @@ -212,6 +303,7 @@ func initGPKGTarget(targetPathFmt string, tmID int, overwrite bool, pagesize int } target := gpkg.TargetGeopackage{} target.Init(targetPath, pagesize) + target.EncodeTiles = encodeTiles return &target } @@ -222,8 +314,109 @@ func injectSuffixIntoPath(p string) string { return path.Join(dir, name+"_%v"+ext) } -func processBySnapping(source processing.Source, targets map[tms20.TMID]processing.Target, tileMatrixSet tms20.TileMatrixSet, snapConfig snap.Config) { - processing.ProcessFeatures(source, targets, func(p geom.Polygon, tmIDs []tms20.TMID) map[tms20.TMID][]geom.Polygon { - return snap.SnapPolygon(p, tileMatrixSet, tmIDs, snapConfig) - }) +func readClip(clipString string) ([]int, error) { + var clip []int + err := json.Unmarshal([]byte(clipString), &clip) + if err != nil { + return nil, err + } + if len(clip) != 3 { + return nil, fmt.Errorf("error: clip must have exactly three values, got %d", len(clip)) + } + return clip, nil +} + +func processBySnapping(source processing.Source, targets map[tms20.TMID]processing.Target, tileMatrixSet tms20.TileMatrixSet, snapConfig processing.Config) { + snapFunction := snap.SnapGeometry + deepestTMID := slices.Max(mapslicehelp.MapKeys(targets)) + rawFactory := pointindex.Factory(tileMatrixSet, deepestTMID) + factory := func() processing.PIndex { return rawFactory() } + processing.ProcessFeatures(source, targets, snapFunction, factory, snapConfig) +} + +// Construct Layer info from config; init gpkg sources +// Make sure to reuse gpkg handles +// Also return helper that closes all sources +func buildLayers(z uint, rawConfig config.TomlConfig) ([]processing.Layer, func()) { + type initedSource struct { + source gpkg.MVTSourceGeopackage + tables map[string]gpkg.Table + } + dataSourceDictionary := processing.DatasourceToDictionary(rawConfig.DataSource) + layers := make([]processing.Layer, 0) + sources := make(map[string]initedSource) + tileset := rawConfig.Tileset + for _, rawLayer := range tileset.Layer { + if z < rawLayer.MinZoom || z > rawLayer.MaxZoom { + continue + } + // Only init gpkg sources once + if _, present := sources[rawLayer.DataSource]; !present { + source, tableMap := initMvtSource(rawLayer, dataSourceDictionary) + sources[rawLayer.DataSource] = initedSource{source, tableMap} + } + initedSource := sources[rawLayer.DataSource] + table, present := initedSource.tables[rawLayer.TableName] + if !present { + log.Printf("table %s not found in datasource %s; skipping layer %s\n", rawLayer.TableName, rawLayer.DataSource, rawLayer.Name) + } else { + // Pair gpkg handle with layer table + source := initedSource.source + source.Table = table + layer := processing.BuildLayer(rawLayer.Name, source) + layers = append(layers, layer) + } + } + closeSources := func() { + for _, s := range sources { + s.source.Close() + } + } + return layers, closeSources +} + +func initMvtSource(rawLayer config.LayerConfig, dataSourceDicationary map[string]string) (gpkg.MVTSourceGeopackage, map[string]gpkg.Table) { + path, present := dataSourceDicationary[rawLayer.DataSource] + if !present { + err := fmt.Errorf("configuration error for layer %s, datasource not found: %s", rawLayer.Name, rawLayer.DataSource) + panic(err) + } + if _, err := os.Stat(path); os.IsNotExist(err) { + err = fmt.Errorf("source Geopackage does not exist: %s", path) + panic(err) + } + source := gpkg.MVTSourceGeopackage{} + source.Init(path) + tables := source.GetTableInfo() + tableMap := make(map[string]gpkg.Table, len(tables)) + for _, table := range tables { + tableMap[table.Name] = table + } + return source, tableMap +} + +// runMVT wires the config, layers and target together to build +// and write the MVT tiles for the requested zoomlevel. +func runMVT(configPath string, zoomlevel uint) error { + if _, err := os.Stat(configPath); os.IsNotExist(err) { + return fmt.Errorf("config file does not exist: %s", configPath) + } + + rawConfig, err := config.ParseMVTConfig(configPath) + if err != nil { + return err + } + if rawConfig.Tileset == nil { + return errors.New("no tileset found in TOML config") + } + + layers, closeSources := buildLayers(zoomlevel, rawConfig) + defer closeSources() + + mvtTarget, err := writetile.NewTileWriter(rawConfig) + if err != nil { + return err + } + + return processing.BuildAndWriteMVTTiles(layers, zoomlevel, mvtTarget) } diff --git a/mapslicehelp/mapslicehelp.go b/mapslicehelp/mapslicehelp.go index 8ce89c1..f06da78 100644 --- a/mapslicehelp/mapslicehelp.go +++ b/mapslicehelp/mapslicehelp.go @@ -51,6 +51,23 @@ func OrderedMapKeys[K comparable, V any](m *orderedmap.OrderedMap[K, V]) []K { return l } +// Assumes all indices 0, 1, ..., len(index)-1 are used exactly once. +func OrderedByIndex[K comparable](index map[K]uint32) []K { + ordered := make([]K, len(index)) + for k, i := range index { + ordered[i] = k + } + return ordered +} + +func MapKeys[K comparable, V any](m map[K]V) []K { + keys := make([]K, 0, len(m)) + for k := range m { + keys = append(keys, k) + } + return keys +} + func RemoveSequences[V any, K comparable](s []V, sequencesToRemove *sortedmap.SortedMap[K, [2]int]) (newS []V) { mmap := sequencesToRemove.Map() keepFrom := 0 diff --git a/pointindex/detect.go b/pointindex/detect.go new file mode 100644 index 0000000..9b5e4dc --- /dev/null +++ b/pointindex/detect.go @@ -0,0 +1,527 @@ +package pointindex + +import ( + "math" + + "github.com/go-spatial/geom" + "github.com/pdok/texel/intgeom" + "github.com/pdok/texel/mapslicehelp" + "github.com/pdok/texel/mathhelp" + "github.com/pdok/texel/morton" + "github.com/pdok/texel/tile" + "github.com/pdok/texel/tms20" +) + +type SegmentIdx struct { + ringIdx int + pointIdx int +} + +type TileClassification int + +const ( + ClassificationUnknown TileClassification = iota + ClassificationIntersect + ClassificationInside + ClassificationOutside +) + +////////////////////////////// +// Tile detection functions // +////////////////////////////// + +func (ix *PointIndex) DetectTilesViaLineTrace(g geom.Geometry, tmsID tms20.TMID, buffer uint) []tile.Tile { + switch g := g.(type) { + case geom.Polygon: + return ix.lineTracePolygon(g, tmsID, buffer) + case geom.LineString: + return ix.lineTraceLine(g, tmsID, buffer) + case geom.Point: + return ix.GetQBBoxWithBuffer(tmsID, buffer) + default: + return nil + + } +} + +// Find all tiles in bbox of geometry currently inserted in index, inflated by buffer +// Goes to pixel level to inflate, then lists all tiles +func (ix *PointIndex) GetQBBoxWithBuffer(tmsID tms20.TMID, bufferSize uint) []tile.Tile { + internalPixelLevel := ix.InternalPixelLevelFromTmsID(tmsID) + tileLevel := levelFromTmsId(tmsID) + quadrants := ix.quadrants[ix.deepestLevel] + + minX := ^uint(0) + maxX := uint(0) + minY := ^uint(0) + maxY := uint(0) + + for z := range quadrants { + x, y := morton.FromZ(z) + minX = min(x, minX) + maxX = max(x, maxX) + minY = min(y, minY) + maxY = max(y, maxY) + } + + var tileMinX, tileMinY uint + if minX < bufferSize { + tileMinX = 0 + } else { + tileMinX = (minX - bufferSize) >> (internalPixelLevel - tileLevel) + } + if minY < bufferSize { + tileMinY = 0 + } else { + tileMinY = (minY - bufferSize) >> (internalPixelLevel - tileLevel) + } + + maxTileCoord := mathhelp.Pow2(tileLevel) - 1 + tileMaxX := min((maxX+bufferSize)>>(internalPixelLevel-tileLevel), maxTileCoord) + tileMaxY := min((maxY+bufferSize)>>(internalPixelLevel-tileLevel), maxTileCoord) + + tiles := make([]tile.Tile, 0, (tileMaxX-tileMinX+1)*(tileMaxY-tileMinY+1)) + for i := range tileMaxX - tileMinX + 1 { + for j := range tileMaxY - tileMinY + 1 { + tileX := tileMinX + i + tileY := tileMinY + j + tiles = append(tiles, ix.makeTile(tileX, tileY, tileLevel, false)) + } + } + return tiles +} + +// Get BBox at deepest level. Currently unused, but more efficient than +// BBox function above when buffer is zero: might optimise in the future. +func (ix *PointIndex) GetPrimitiveQBBox(l Level) []tile.Tile { + quadrants := ix.quadrants[l] + + minX := ^uint(0) + maxX := uint(0) + minY := ^uint(0) + maxY := uint(0) + + for z := range quadrants { + x, y := morton.FromZ(z) + minX = min(x, minX) + maxX = max(x, maxX) + minY = min(y, minY) + maxY = max(y, maxY) + } + + quadrantSlice := make([]tile.Tile, (maxY-minY+1)*(maxX-minX+1)) + for i := range maxX - minX + 1 { + for j := range maxY - minY + 1 { + extent, _ := ix.getQuadrantExtentAndCentroid(l, minX+i, minY+j, ix.intExtent) + newTile := tile.Tile{ + Extent: extent, + X: minX + i, + Y: minY + i, + IsContained: false, + } + quadrantSlice[i*(maxY-minY+1)+j] = newTile + } + } + return quadrantSlice +} + +////////////////////////// +// General line tracing // +////////////////////////// + +// registerFunc marks a tile at (xCoord, yCoord) (in tile coordinates at +// level l) as touched. Callers that need to associate the touched tile +// with a segment should do so via a variable captured in their closure. +type registerFunc func(xCoord, yCoord uint, l Level) + +// Line trace over a line, registering tiles hit with `register`. +// Given a line, detect in which quadrant the line advances. Essentially do a breadth-first search in that direction. +func (ix *PointIndex) lineTrace(line geom.Line, tileLevel, intPixLevel Level, buffer uint, register registerFunc) { + intLine := intgeom.FromGeomLine(line) + + // Detect direction + var dx, dy int + if intLine.Point1().X() < intLine.Point2().X() { + dx = 1 + } else { + dx = -1 + } + if intLine.Point1().Y() < intLine.Point2().Y() { + dy = 1 + } else { + dy = -1 + } + + startTileX, startTileY := ix.findTile(intLine.Point1(), tileLevel) + startX, startY := int(startTileX), int(startTileY) //nolint:gosec // G115 + bufferSize := ix.getResolution(intPixLevel) * intgeom.M(buffer) //nolint:gosec // G115 + + // We will advance in a certain direction. + // Here, we cover the tiles "behind" us. + ix.tryRegisterTile(intLine, startX-dx, startY+dy, tileLevel, bufferSize, register) + ix.tryRegisterTile(intLine, startX-dx, startY, tileLevel, bufferSize, register) + ix.tryRegisterTile(intLine, startX-dx, startY-dy, tileLevel, bufferSize, register) + ix.tryRegisterTile(intLine, startX, startY-dy, tileLevel, bufferSize, register) + ix.tryRegisterTile(intLine, startX+dx, startY-dy, tileLevel, bufferSize, register) + ix.tryRegisterTile(intLine, startX, startY, tileLevel, bufferSize, register) + + type coord struct{ x, y int } + + // Advance in steps. At each step, check all relevant tiles of an anti-diagonal. + // Initial anti-diagonal contains only start tile. + // We keep track of the previous two steps; we need one extra step for the rare case a + // line passes exactly through a corner at buffer = 0. + frontier := []coord{{startX, startY}} + prevFrontier := make([]coord, 0) + + // Candidates for the next step are determined by previous successes. Loop stops after we have + // run out of successes. + for len(prevFrontier)+len(frontier) > 0 { + // Set of next candidates + nextSet := make(map[coord]bool, len(frontier)*2+len(prevFrontier)) + for _, cur := range frontier { + nextSet[coord{cur.x + dx, cur.y}] = true + nextSet[coord{cur.x, cur.y + dy}] = true + } + for _, prev := range prevFrontier { + nextSet[coord{prev.x + dx, prev.y + dy}] = true + } + + // Try registering and record successes + next := make([]coord, 0, len(nextSet)) + for c := range nextSet { + if ix.tryRegisterTile(intLine, c.x, c.y, tileLevel, bufferSize, register) { + next = append(next, c) + } + } + + // Update set of past successes + prevFrontier = frontier + frontier = next + } +} + +// Check whether tile hits line and registers upon success. Also returns success. +func (ix *PointIndex) tryRegisterTile(line intgeom.Line, x, y int, l Level, bufferSize intgeom.M, register registerFunc) bool { + maxCoord := int(mathhelp.Pow2(l)) - 1 //nolint:gosec // G115 + if x < 0 || y < 0 || x > maxCoord || y > maxCoord { + return false // out of bounds: no tile there, nothing to register + } + ux, uy := uint(x), uint(y) + extent, _ := ix.getQuadrantExtentAndCentroid(l, ux, uy, ix.intExtent) + if tileIntersectsLine(line, extent, bufferSize) { + register(ux, uy, l) + return true + } + return false +} + +// Detect whether line intersects extent, inflated by buffer. +// Uses same logic as quadrants: left and bottom edge belong to extent, top and right do not. +func tileIntersectsLine(line intgeom.Line, extent intgeom.Extent, buffer intgeom.M) bool { + bufferedExtent := intgeom.Extent{ + extent.MinX() - buffer, + extent.MinY() - buffer, + extent.MaxX() + buffer, + extent.MaxY() + buffer, + } + + return lineIntersects(line, bufferedExtent) +} + +/////////////////////////////// +// Line tracing for polygons // +/////////////////////////////// + +// Line tracing for polygons is special in the sense that we do not just trace the lines, +// but also find tiles that lie entirely on the interior of the polygon. This second process +// we call "filling". + +// Detect tiles for polygon using line tracing. +// First line trace, then fill remaining tiles. This enables us to classify tiles as being "Inside" and "Outside" the polygon. +// Then process output of `fillTiles` to return a list at tile level. +func (ix *PointIndex) lineTracePolygon(polygon geom.Polygon, tmsID tms20.TMID, buffer uint) (tiles []tile.Tile) { + tileLevel := levelFromTmsId(tmsID) + segements, classification := ix.lineTracePolygonEdges(polygon, tmsID, buffer) + ix.fillTiles(tileLevel, 0, 0, true, segements, classification, polygon) + + // Extract list of tiles at tile level. + tiles = make([]tile.Tile, 0) + for level, classificationAtLevel := range classification { + for z, class := range classificationAtLevel { + switch class { + case ClassificationInside: + size := mathhelp.Pow2((tileLevel - level) * 2) + baseZ := z << (2 * (tileLevel - level)) + for i := range size { + tiles = append(tiles, ix.makeTileZ(baseZ+i, tileLevel, true)) + } + case ClassificationOutside: + continue + case ClassificationIntersect: + if level == tileLevel { + tiles = append(tiles, ix.makeTileZ(z, tileLevel, false)) + } + case ClassificationUnknown: + panic("ClassificationUnknown tile for polygon during linetrace") + } + } + } + return tiles +} + +// Loop over polygon edges and apply line trace. +// Use a special registering function that keeps track of which edges hit which tiles. +func (ix *PointIndex) lineTracePolygonEdges(polygon geom.Polygon, tmsID tms20.TMID, buffer uint) (segments map[morton.Z][]SegmentIdx, classification map[Level]map[morton.Z]TileClassification) { + tileLevel := levelFromTmsId(tmsID) + intPixLevel := ix.InternalPixelLevelFromTmsID(tmsID) + + // Initialize data + segments = make(map[morton.Z][]SegmentIdx) + classification = make(map[Level]map[morton.Z]TileClassification, tileLevel+1) + for l := range tileLevel + 1 { + classification[l] = make(map[morton.Z]TileClassification) + } + + // Helper registering function that marks a tile and all its parents as "Intersecting". + var markIntersected func(l Level, z morton.Z) + markIntersected = func(l Level, z morton.Z) { + if classification[l][z] == ClassificationIntersect { + return + } + classification[l][z] = ClassificationIntersect + if l == 0 { + return + } + markIntersected(l-1, z>>2) + } + + // Register function that registers when a segment hits a tile and then marks tiles as "Intersected". + // Note the use of `idx` in this closure: we manually update it in the loop below. + var idx SegmentIdx + register := func(x, y uint, l Level) { + z := morton.MustToZ(x, y) + segments[z] = append(segments[z], idx) + markIntersected(l, z) + } + + // Loop over all edges + for ringIdx, ring := range polygon.LinearRings() { + for pointIdx := range ring { + // Update `idx` for use in the `register` closure + idx = SegmentIdx{ringIdx: ringIdx, pointIdx: pointIdx} + // Apply line trace + line := geom.Line{ring[pointIdx], ring[(pointIdx+1)%len(ring)]} + ix.lineTrace(line, tileLevel, intPixLevel, buffer, register) + } + } + return segments, classification +} + +// Recurse into quadtree, filling entire area as either "Inside", "Outside", or "Intersect" at tile level. +// A classification of "Inside" or "Outside" at a higher level means that all children at tile level have this classification. +// We assume that "Intersect" at the tile level is already classified, and that parents are also marked as "Intersect". +// This means that tiles marked as "Intersect" will need to be recursed into for fine-grained classification. +// It also means that at any level, all children of an unmarked tile have the same classification +// This enables the classification at higher levels. +// We also assume that `segments` contains the data of which segments of polygon intersect which tiles at tile level. +func (ix *PointIndex) fillTiles(tileLevel, currentLevel Level, currentZ morton.Z, containsAll bool, segments map[morton.Z][]SegmentIdx, classification map[Level]map[morton.Z]TileClassification, polygon geom.Polygon) { + if tileLevel == currentLevel { + return + } + + // Get children + children := getChildren(currentZ) + intersectingChildren := make([]morton.Z, 0, 4) + nextLevel := currentLevel + 1 + + // First classify unknown children. This can be done at a higher level + for _, child := range children { + if _, present := classification[nextLevel][child]; present { + intersectingChildren = append(intersectingChildren, child) + continue + } + if containsAll { + classification[nextLevel][child] = ClassificationOutside + } else { + childAtTileLevel := child << ((tileLevel - nextLevel) * 2) + classification[nextLevel][child] = ix.fillTile(childAtTileLevel, tileLevel, segments, classification, polygon) + } + } + + containsAll = containsAll && len(intersectingChildren) < 2 + + // For intersecting children we recurse to the next level + for _, child := range intersectingChildren { + if _, present := classification[nextLevel][child]; present { + ix.fillTiles(tileLevel, nextLevel, child, containsAll, segments, classification, polygon) + } + } +} + +func getChildren(z morton.Z) [4]morton.Z { + shift := z << 2 + return [4]morton.Z{shift, shift + 1, shift + 2, shift + 3} +} + +// Use raycast to classify a tile as "Inside" or "Outside". Does not work on tiles that are "Intersect". +// tileLevel needs to be the actual tile level and z needs to be a coordinate at this level +// Relies on classification and segments as pre-computed values. +func (ix *PointIndex) fillTile(z morton.Z, tileLevel Level, segments map[morton.Z][]SegmentIdx, classification map[Level]map[morton.Z]TileClassification, polygon geom.Polygon) TileClassification { + x, y := morton.FromZ(z) + intersectingTilesLeft := ix.findIntersectingTilesLeft(x, y, tileLevel, classification) + + tileHeightCoord := ix.getResolution(tileLevel)*intgeom.M(y) + ix.intExtent.MinY() //nolint:gosec // G115 + + seen := make(map[SegmentIdx]bool) + numIntersections := 0 + + // Raycast and loop over intersecting segments + for _, z := range intersectingTilesLeft { + for _, segment := range segments[z] { + if seen[segment] { + continue + } + seen[segment] = true + ring := polygon.LinearRings()[segment.ringIdx] + + y1 := intgeom.FromGeomOrd(ring[segment.pointIdx][1]) + y2 := intgeom.FromGeomOrd(ring[(segment.pointIdx+1)%len(ring)][1]) + + minY := min(y1, y2) + maxY := max(y1, y2) + + // Logic for deciding relevant segments + // Needed for when ray intersects vertex of polygon + switch { + case minY == maxY: + case maxY < tileHeightCoord: + case minY >= tileHeightCoord: + default: + numIntersections++ + } + } + } + + // Result determined by "Jordan Curve Theorem" + if numIntersections%2 == 0 { + return ClassificationOutside + } + return ClassificationInside +} + +// Helper function for raycasting function `fillTile`. Finds all tiles left of argument containing segments of polygon. +// tileLevel needs to be the actual tile level and x and y need to be coordinates for this level. +// Relies on `classification` as pre-computed information. +// Searches by binary search. +func (ix *PointIndex) findIntersectingTilesLeft(x, y, tileLevel Level, classification map[Level]map[morton.Z]TileClassification) []morton.Z { + intersectingCurrentLevel := []morton.Z{0} + var intersectingNextLevel []morton.Z + var leftChild, rightChild morton.Z + + // Binary search over the tile level + for currentLevel := range tileLevel { + intersectingNextLevel = make([]morton.Z, 0) + + // Coordinates of parent of tile at the next level + xAtNextLevel := x >> (tileLevel - currentLevel - 1) + yAtNextLevel := y >> (tileLevel - currentLevel - 1) + + // Decide whether we need the upper or lower children for the next level + nextLevelDown := yAtNextLevel%2 == 0 + for _, z := range intersectingCurrentLevel { + // Compute candidates + if nextLevelDown { + leftChild = z << 2 + rightChild = (z << 2) + 1 + } else { + leftChild = (z << 2) + 2 + rightChild = (z << 2) + 3 + } + + // If left tile present, register it + if _, present := classification[currentLevel+1][leftChild]; present { + intersectingNextLevel = append(intersectingNextLevel, leftChild) + } + + // Right tile needs to be present AND to the left of argument + rightX, _ := morton.FromZ(rightChild) + if _, present := classification[currentLevel+1][rightChild]; present && rightX <= xAtNextLevel { + intersectingNextLevel = append(intersectingNextLevel, rightChild) + } + } + intersectingCurrentLevel = intersectingNextLevel + } + return intersectingCurrentLevel +} + +//////////////////////////// +// Line tracing for lines // +//////////////////////////// + +// Line trace by looping over segments +func (ix *PointIndex) lineTraceLine(line geom.LineString, tmsID tms20.TMID, buffer uint) []tile.Tile { + tileSet := make(map[tile.Tile]bool) + tileLevel := levelFromTmsId(tmsID) + internalPixelLevel := ix.InternalPixelLevelFromTmsID(tmsID) + + register := func(x, y uint, l Level) { + tile := ix.makeTile(x, y, l, false) + tileSet[tile] = true + } + + segments, _ := line.AsSegments() + for _, segment := range segments { + ix.lineTrace(segment, tileLevel, internalPixelLevel, buffer, register) + } + + return mapslicehelp.MapKeys(tileSet) +} + +////////////////////// +// Helper functions // +////////////////////// + +func (ix *PointIndex) makeTileZ(z morton.Z, l Level, isContained bool) tile.Tile { + x, y := morton.FromZ(z) + return ix.makeTile(x, y, l, isContained) +} + +func (ix *PointIndex) makeTile(x, y uint, l Level, isContained bool) tile.Tile { + extent, _ := ix.getQuadrantExtentAndCentroid(l, x, y, ix.intExtent) + return tile.Tile{ + Extent: extent, + X: x, + Y: mathhelp.Pow2(l) - 1 - y, + IsContained: isContained, + } +} + +func (ix *PointIndex) getResolution(level Level) intgeom.M { + return ix.intExtent.XSpan() / int64(mathhelp.Pow2(level)) //nolint:gosec // G115 +} + +func (ix *PointIndex) InternalPixelLevelFromTmsID(deepestTIMID tms20.TMID) Level { + levelDiff := uint(math.Log2(float64(ix.tilePixels))) + uint(math.Log2(float64(ix.internalPixels))) + return uint(deepestTIMID) + levelDiff //nolint:gosec // G115 +} + +func (ix *PointIndex) TmsIDFromInternalPixelLevel(level Level) tms20.TMID { + levelDiff := uint(math.Log2(float64(ix.tilePixels))) + uint(math.Log2(float64(ix.internalPixels))) + return int(level) - int(levelDiff) //nolint:gosec // G115 +} + +func levelFromTmsId(tms tms20.TMID) Level { + return Level(tms) //nolint:gosec // G115 integers < 40 +} + +func (ix *PointIndex) findTile(p *intgeom.Point, l Level) (tileX, tileY uint) { + levelDiff := ix.deepestLevel - l + //nolint:gosec // G115 + deepestTileX := uint(p.X()-ix.intExtent.MinX()) / uint(ix.deepestRes) + //nolint:gosec // G115 + deepestTileY := uint(p.Y()-ix.intExtent.MinY()) / uint(ix.deepestRes) + tileX = deepestTileX >> levelDiff + tileY = deepestTileY >> levelDiff + return +} diff --git a/pointindex/detect_test.go b/pointindex/detect_test.go new file mode 100644 index 0000000..eb8481c --- /dev/null +++ b/pointindex/detect_test.go @@ -0,0 +1,789 @@ +package pointindex + +import ( + "sort" + "testing" + + "github.com/go-spatial/geom" + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" + + "github.com/pdok/texel/intgeom" + "github.com/pdok/texel/mathhelp" + "github.com/pdok/texel/morton" + "github.com/pdok/texel/tile" + "github.com/pdok/texel/tms20" +) + +// Note that wantTiles is specified with left-down origin, but makeTile introduces the left-up origin. +func TestPointIndex_GetQBBoxWithBuffer(t *testing.T) { + type tileCoord struct{ x, y uint } + tests := []struct { + name string + tmsID tms20.TMID + tilePixels uint + deepestLevel Level + points [][2]int // Points to be inserted at deepestlevel + bufferSize uint + wantTiles []tileCoord + }{ + { + // deepestLevel 4 -> 16x16 pixel grid, tmsID 2 -> 4x4 tiles of 4x4 pixels each. + name: "points all within one tile, no buffer", + tmsID: 2, + tilePixels: 4, + points: [][2]int{{5, 5}, {6, 6}}, + bufferSize: 0, + wantTiles: []tileCoord{{1, 1}}, + }, + { + name: "point in one tile, buffer expands to multiple tiles", + tmsID: 2, + tilePixels: 4, + points: [][2]int{{0, 0}}, + bufferSize: 6, + wantTiles: []tileCoord{{0, 0}, {0, 1}, {1, 0}, {1, 1}}, + }, + { + name: "non-square rectangle of tiles", + tmsID: 2, + tilePixels: 4, + points: [][2]int{{4, 4}, {11, 5}}, + bufferSize: 0, + wantTiles: []tileCoord{{1, 1}, {2, 1}}, + }, + { + name: "points in diagonal tiles, expands to rectangle", + tmsID: 2, + tilePixels: 4, + points: [][2]int{{3, 3}, {4, 4}}, + bufferSize: 0, + wantTiles: []tileCoord{{0, 0}, {0, 1}, {1, 0}, {1, 1}}, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + ix := newPointIndex(tt.tmsID, tt.tilePixels, 1, 1.0) + for _, p := range tt.points { + require.NoError(t, ix.InsertCoord(p[0], p[1])) + } + + l := levelFromTmsId(tt.tmsID) + want := make([]tile.Tile, 0, len(tt.wantTiles)) + for _, tc := range tt.wantTiles { + want = append(want, ix.makeTile(tc.x, tc.y, l, false)) + } + + got := ix.GetQBBoxWithBuffer(tt.tmsID, tt.bufferSize) + assert.ElementsMatch(t, want, got) + }) + } +} + +// registeredTile records a single call to a RegisterFunc during a test. +type registeredTile struct { + x, y uint +} + +// Register function for testing +func recordingRegister(dst *[]registeredTile) registerFunc { + return func(tileX, tileY uint, l Level) { + maxCoord := mathhelp.Pow2(l) - 1 + if tileX > maxCoord || tileY > maxCoord { + return + } + *dst = append(*dst, registeredTile{tileX, tileY}) + } +} + +// order and deduplicate slice of registeredTile (for testing) +func uniqueTileCoords(records []registeredTile) [][2]uint { + seen := map[[2]uint]bool{} + var out [][2]uint + for _, r := range records { + key := [2]uint{r.x, r.y} + if !seen[key] { + seen[key] = true + out = append(out, key) + } + } + sort.Slice(out, func(i, j int) bool { + if out[i][0] != out[j][0] { + return out[i][0] < out[j][0] + } + return out[i][1] < out[j][1] + }) + return out +} + +// Pointindex with centre not at (0,0) +func newOffsetPointIndex(deepestLevel Level, cellSize, originX, originY float64) *PointIndex { + deepestSize := mathhelp.Pow2(deepestLevel) + span := cellSize * float64(deepestSize) + intExtent := intgeom.Extent{ + intgeom.FromGeomOrd(originX), intgeom.FromGeomOrd(originY), + intgeom.FromGeomOrd(originX + span), intgeom.FromGeomOrd(originY + span), + } + return &PointIndex{ + Quadrant: Quadrant{intExtent: intExtent}, + deepestLevel: deepestLevel, + deepestSize: deepestSize, + //nolint:gosec // G115 + deepestRes: intExtent.XSpan() / int64(deepestSize), + } +} + +func TestLineTrace_TilesTouched(t *testing.T) { + tests := []struct { + name string + deepestLevel Level + l Level + cellSize float64 + originX float64 + originY float64 + line geom.Line + buffer uint + want [][2]uint + }{ + { + name: "no buffer, +x+y diagonal", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 2}, geom.Point{14, 10}}, + buffer: 0, + want: [][2]uint{{0, 0}, {1, 0}, {1, 1}, {2, 1}, {2, 2}, {3, 2}}, + }, + { + name: "no buffer, -x-y diagonal (reverse of +x+y)", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{14, 10}, geom.Point{2, 2}}, + buffer: 0, + want: [][2]uint{{0, 0}, {1, 0}, {1, 1}, {2, 1}, {2, 2}, {3, 2}}, + }, + { + name: "no buffer, +x-y diagonal", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 10}, geom.Point{14, 2}}, + buffer: 0, + want: [][2]uint{{0, 2}, {1, 1}, {1, 2}, {2, 0}, {2, 1}, {3, 0}}, + }, + { + name: "no buffer, -x+y diagonal (reverse of +x-y)", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{14, 2}, geom.Point{2, 10}}, + buffer: 0, + want: [][2]uint{{0, 2}, {1, 1}, {1, 2}, {2, 0}, {2, 1}, {3, 0}}, + }, + { + name: "buffer, +x+y diagonal (reaches beyond index limit)", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 2}, geom.Point{14, 9}}, + buffer: 2, + want: [][2]uint{{0, 0}, {0, 1}, {1, 0}, {1, 1}, {1, 2}, {2, 0}, {2, 1}, {2, 2}, {3, 1}, {3, 2}}, + }, + { + name: "buffer, -x-y diagonal (ververse of +x+y)", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{14, 9}, geom.Point{2, 2}}, + buffer: 2, + want: [][2]uint{{0, 0}, {0, 1}, {1, 0}, {1, 1}, {1, 2}, {2, 0}, {2, 1}, {2, 2}, {3, 1}, {3, 2}}, + }, + { + name: "buffer, +x-y diagonal", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 10}, geom.Point{14, 2}}, + buffer: 2, + want: [][2]uint{{0, 1}, {0, 2}, {0, 3}, {1, 0}, {1, 1}, {1, 2}, {1, 3}, {2, 0}, {2, 1}, {2, 2}, {3, 0}, {3, 1}}, + }, + { + name: "buffer, -x+y diagonal: (reverse of +x-y)", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{14, 2}, geom.Point{2, 10}}, + buffer: 2, + want: [][2]uint{{0, 1}, {0, 2}, {0, 3}, {1, 0}, {1, 1}, {1, 2}, {1, 3}, {2, 0}, {2, 1}, {2, 2}, {3, 0}, {3, 1}}, + }, + { + name: "horizontal line", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 5}, geom.Point{14, 5}}, + buffer: 0, + want: [][2]uint{{0, 1}, {1, 1}, {2, 1}, {3, 1}}, + }, + { + name: "vertical line", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{5, 2}, geom.Point{5, 14}}, + buffer: 0, + want: [][2]uint{{1, 0}, {1, 1}, {1, 2}, {1, 3}}, + }, + { + name: "end in upper right corner", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 2}, geom.Point{3, 3}}, + buffer: 0, + want: [][2]uint{{0, 0}}, + }, + { + name: "alongside upper edge", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 2}, geom.Point{3, 3}}, + buffer: 1, + want: [][2]uint{{0, 0}, {0, 1}, {1, 0}, {1, 1}}, + }, + { + name: "start in corner, register tiles behind you", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{3, 3}, geom.Point{2, 2}}, + buffer: 1, + want: [][2]uint{{0, 0}, {0, 1}, {1, 0}, {1, 1}}, + }, + { + name: "no buffer: walk across tile corner", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{6, 6}, geom.Point{10, 10}}, + buffer: 0, + want: [][2]uint{{1, 1}, {1, 2}, {2, 1}, {2, 2}}, + }, + { + name: "buffer: walk across tile corner", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{6, 6}, geom.Point{7, 7}}, + buffer: 1, + want: [][2]uint{{1, 1}, {1, 2}, {2, 1}, {2, 2}}, + }, + { + name: "buffer clips upper right corner of tile", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{12, 14}, geom.Point{14, 12}}, + buffer: 1, + want: [][2]uint{{2, 2}, {2, 3}, {3, 2}, {3, 3}}, // Arguably {2,2} should not be here + }, + { + name: "non-zero origin: same relative result as +x+y diagonal", + deepestLevel: 4, l: 2, cellSize: 1.0, originX: 100, originY: 100, + line: geom.Line{geom.Point{102, 102}, geom.Point{114, 110}}, + buffer: 0, + want: [][2]uint{{0, 0}, {1, 0}, {1, 1}, {2, 1}, {2, 2}, {3, 2}}, + }, + { + name: "non-zero origin with buffer: same relative result as buffered corner case", + deepestLevel: 4, l: 2, cellSize: 1.0, originX: 100, originY: 100, + line: geom.Line{geom.Point{103, 103}, geom.Point{105, 105}}, + buffer: 1, + want: [][2]uint{{0, 0}, {0, 1}, {1, 0}, {1, 1}}, + }, + } + + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + ix := newOffsetPointIndex(tt.deepestLevel, tt.cellSize, tt.originX, tt.originY) + + var recorded []registeredTile + ix.lineTrace(tt.line, tt.l, tt.deepestLevel, tt.buffer, recordingRegister(&recorded)) + + got := uniqueTileCoords(recorded) + assert.Equal(t, tt.want, got) + }) + } +} + +// helper for building classification data +func buildClassification(tuples ...[3]uint) map[Level]map[morton.Z]TileClassification { + classification := make(map[Level]map[morton.Z]TileClassification) + for _, tuple := range tuples { + level, x, y := tuple[0], tuple[1], tuple[2] + if classification[level] == nil { + classification[level] = make(map[morton.Z]TileClassification) + } + classification[level][morton.MustToZ(x, y)] = ClassificationIntersect + } + return classification +} + +func TestPointIndex_findIntersectingTilesLeft(t *testing.T) { + tests := []struct { + name string + x, y uint + targetLevel Level + classification map[Level]map[morton.Z]TileClassification + want []morton.Z + }{ + { + name: "trivial example at level 0", + x: 0, + y: 0, + targetLevel: 0, + classification: buildClassification(), + want: []morton.Z{0}, + }, + { + name: "no tiles", + x: 3, + y: 3, + targetLevel: 2, + classification: buildClassification(), + want: []morton.Z{}, + }, + { + name: "single tile registered, check left from this tile.", + x: 0, + y: 0, + targetLevel: 3, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 0, 0}, + [3]uint{2, 0, 0}, + [3]uint{3, 0, 0}, + ), + want: []morton.Z{0}, + }, + { + name: "single tile registered, check this tile, requires rightChild", + x: 1, + y: 0, + targetLevel: 1, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 1, 0}, + ), + want: []morton.Z{1}, + }, + { + name: "tile registered right of target", + x: 0, + y: 0, + targetLevel: 1, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 1, 0}, + ), + want: []morton.Z{}, + }, + { + name: "use odd y coordinate", + x: 0, + y: 1, + targetLevel: 1, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 0, 1}, + ), + want: []morton.Z{2}, + }, + { + name: "use odd y and rightChild", + x: 0, + y: 1, + targetLevel: 1, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 1, 1}, + ), + want: []morton.Z{}, + }, + { + name: "ordd y and rightChild, want to find it", + x: 1, + y: 1, + targetLevel: 1, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 1, 1}, + ), + want: []morton.Z{3}, + }, + { + name: "register entire row", + x: 3, + y: 0, + targetLevel: 2, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 0, 0}, + [3]uint{1, 1, 0}, + [3]uint{2, 0, 0}, + [3]uint{2, 1, 0}, + [3]uint{2, 2, 0}, + [3]uint{2, 3, 0}, + ), + want: []morton.Z{0, 1, 4, 5}, + }, + { + name: "corners of grid populated", + x: 7, + y: 7, + targetLevel: 3, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 0, 0}, [3]uint{1, 0, 1}, [3]uint{1, 1, 0}, [3]uint{1, 1, 1}, + [3]uint{2, 0, 0}, [3]uint{2, 0, 3}, [3]uint{2, 3, 0}, [3]uint{2, 3, 3}, + [3]uint{3, 0, 0}, [3]uint{3, 7, 0}, [3]uint{3, 0, 7}, [3]uint{3, 7, 7}, + ), + want: []morton.Z{42, 63}, + }, + { + name: "general test with two disjoint tiles to pick up", + x: 6, + y: 7, + targetLevel: 3, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 0, 1}, [3]uint{1, 1, 1}, + [3]uint{2, 0, 3}, [3]uint{2, 2, 3}, [3]uint{2, 3, 3}, + [3]uint{3, 0, 7}, [3]uint{2, 6, 1}, [3]uint{3, 5, 7}, [3]uint{3, 7, 7}, + ), + want: []morton.Z{42, 59}, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + ix := newSimplePointIndex(tt.targetLevel, 1.0) + got := ix.findIntersectingTilesLeft(tt.x, tt.y, tt.targetLevel, tt.classification) + t.Logf("findIntersectingTilesLeft(%d, %d, %d) = %v", tt.x, tt.y, tt.targetLevel, got) + assert.Equal(t, tt.want, got) + }) + } +} + +func squarePolygon(minX, minY, maxX, maxY float64) geom.Polygon { + return geom.Polygon{ + {{minX, minY}, {maxX, minY}, {maxX, maxY}, {minX, maxY}}, + } +} + +func squareWithHolePolygon(minX, minY, maxX, maxY, holeMinX, holeMinY, holeMaxX, holeMaxY float64) geom.Polygon { + return geom.Polygon{ + {{minX, minY}, {maxX, minY}, {maxX, maxY}, {minX, maxY}}, + {{holeMinX, holeMinY}, {holeMinX, holeMaxY}, {holeMaxX, holeMaxY}, {holeMaxX, holeMinY}}, + } +} + +func trianglePolygon(x1, y1, x2, y2, x3, y3 float64) geom.Polygon { + return geom.Polygon{ + {{x1, y1}, {x2, y2}, {x3, y3}}, + } +} + +func TestPointIndex_classifyNonIntersectingTile(t *testing.T) { + tests := []struct { + name string + ix *PointIndex + polygon geom.Polygon + tmsID tms20.TMID + buffer uint + tileLevel Level + tileX, tileY uint + want TileClassification + }{ + { + name: "tile outside a centered square", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(2, 2, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 0, + tileY: 0, + want: ClassificationOutside, + }, + { + name: "tile outside a centered square", + ix: newSimplePointIndex(3, 1.5), + polygon: squarePolygon(3, 3, 9, 9), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 0, + tileY: 0, + want: ClassificationOutside, + }, + { + name: "tile inside a centered square", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(2, 2, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 3, + tileY: 3, + want: ClassificationInside, + }, + { + name: "higher-level tile inside a centered square", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(1, 1, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 2, + tileX: 1, + tileY: 1, + want: ClassificationInside, + }, + { + name: "tile clearly outside, opposite corner", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(2, 2, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 7, + tileY: 7, + want: ClassificationOutside, + }, + { + name: "tile inside the hole of a donut polygon", + ix: newSimplePointIndex(3, 1.0), + // This should work with hole 3,3,5,5 but lineIntesects is buggy + polygon: squareWithHolePolygon(0, 0, 7, 7, 3, 3, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 4, + tileY: 4, + want: ClassificationOutside, + }, + { + name: "tile in the solid part of a donut polygon", + ix: newSimplePointIndex(3, 1.0), + polygon: squareWithHolePolygon(0, 0, 8, 8, 3, 3, 5, 5), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 1, + tileY: 1, + want: ClassificationInside, + }, + { + name: "tile in the solid part of a donut polygon, far corner", + ix: newSimplePointIndex(3, 1.0), + polygon: squareWithHolePolygon(0, 0, 8, 8, 3, 3, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 7, + tileY: 7, + want: ClassificationOutside, + }, + { + name: "tile below a diagonal triangle (raycast hits intersection)", + ix: newSimplePointIndex(3, 1.0), + polygon: trianglePolygon(0, 0, 0, 7, 7, 7), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 1, + tileY: 0, + want: ClassificationOutside, + }, + { + name: "tile above a diagonal triangle (raycast hits intersection)", + ix: newSimplePointIndex(3, 1.0), + polygon: trianglePolygon(0, 0, 7, 0, 0, 7), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 1, + tileY: 7, + want: ClassificationOutside, + }, + { + name: "nonzero buffer around a centered square", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(2, 2, 6, 6), + tmsID: 3, + buffer: 4, + tileLevel: 3, + tileX: 0, + tileY: 0, + want: ClassificationOutside, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + segments, classification := tt.ix.lineTracePolygonEdges(tt.polygon, tt.tmsID, tt.buffer) + z := morton.MustToZ(tt.tileX, tt.tileY) + got := tt.ix.fillTile(z, tt.tileLevel, segments, classification, tt.polygon) + t.Logf("classifyNonIntersectingTile(tile=(%d,%d)) = %v", tt.tileX, tt.tileY, got) + assert.Equal(t, tt.want, got) + }) + } +} + +// convert classification data to visual test data +func classificationGrid(classification map[Level]map[morton.Z]TileClassification, targetLevel Level) map[Level][][]TileClassification { + grid := make(map[Level][][]TileClassification, targetLevel+1) + for l := Level(0); l <= targetLevel; l++ { + levelClassification := classification[l] + size := uint(1) << l + rows := make([][]TileClassification, size) + for y := range size { + row := make([]TileClassification, size) + for x := range size { + z := morton.MustToZ(x, y) + if c, ok := levelClassification[z]; ok { + row[x] = c + } else { + row[x] = ClassificationUnknown + } + } + rows[y] = row + } + grid[l] = rows + } + return grid +} + +func TestPointIndex_classifyNonIntersectingTiles(t *testing.T) { + const ( + u = ClassificationUnknown + x = ClassificationIntersect + i = ClassificationInside + o = ClassificationOutside + ) + tests := []struct { + name string + ix *PointIndex + polygon geom.Polygon + tmsID tms20.TMID + buffer uint + want map[Level][][]TileClassification + }{ + { + name: "centered square", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(2, 2, 6, 6), + tmsID: 3, + buffer: 0, + want: map[Level][][]TileClassification{ + 0: {{x}}, + 1: {{x, x}, {x, x}}, + 2: { + {o, x, x, o}, + {x, x, x, x}, + {x, x, x, x}, + {o, x, x, x}, + }, + 3: { + {u, u, o, o, o, o, u, u}, + {u, u, x, x, x, x, u, u}, + {o, x, x, x, x, x, x, o}, + {o, x, x, i, i, x, x, o}, + {o, x, x, i, i, x, x, o}, + {o, x, x, x, x, x, x, o}, + {u, u, x, x, x, x, x, o}, + {u, u, o, o, o, o, o, o}, + }, + }, + }, + { + name: "donut polygon (square with a square hole)", + ix: newSimplePointIndex(3, 1.0), + polygon: squareWithHolePolygon(0, 0, 7, 7, 3, 3, 5, 5), + tmsID: 3, + buffer: 0, + want: map[Level][][]TileClassification{ + 0: {{x}}, + 1: {{x, x}, {x, x}}, + 2: { + {x, x, x, x}, + {x, x, x, x}, + {x, x, x, x}, + {x, x, x, x}, + }, + 3: { + {x, x, x, x, x, x, x, x}, + {x, i, i, i, i, i, x, x}, + {x, i, i, x, x, i, x, x}, + {x, i, x, x, x, x, x, x}, + {x, i, x, x, x, x, x, x}, + {x, i, i, x, x, x, x, x}, + {x, x, x, x, x, x, x, x}, + {x, x, x, x, x, x, x, x}, + }, + }, + }, + { + name: "diagonal triangle", + ix: newSimplePointIndex(3, 1.0), + polygon: trianglePolygon(0, 0, 7, 0, 7, 7), + tmsID: 3, + buffer: 0, + want: map[Level][][]TileClassification{ + 0: {{x}}, + 1: {{x, x}, {x, x}}, + 2: { + {x, x, x, x}, + {x, x, x, x}, + {o, x, x, x}, + {o, o, x, x}, + }, + 3: { + {x, x, x, x, x, x, x, x}, + {x, x, x, i, i, i, x, x}, + {o, x, x, x, i, i, x, x}, + {o, o, x, x, x, i, x, x}, + {u, u, o, x, x, x, x, x}, + {u, u, o, o, x, x, x, x}, + {u, u, u, u, o, x, x, x}, + {u, u, u, u, o, o, o, x}, + }, + }, + }, + { + name: "polygon covering the whole extent", + ix: newSimplePointIndex(2, 1.0), + polygon: squarePolygon(0, 0, 4, 4), + tmsID: 2, + buffer: 0, + want: map[Level][][]TileClassification{ + 0: {{x}}, + 1: {{x, x}, {x, x}}, + 2: { + {x, x, x, x}, + {x, i, i, x}, + {x, i, i, x}, + {x, x, x, x}, + }, + }, + }, + { + name: "tiny polygon in a single corner tile", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(0.25, 0.25, 0.75, 0.75), + tmsID: 3, + buffer: 0, + want: map[Level][][]TileClassification{ + 0: {{x}}, + 1: {{x, o}, {o, o}}, + 2: { + {x, o, u, u}, + {o, o, u, u}, + {u, u, u, u}, + {u, u, u, u}, + }, + 3: { + {x, o, u, u, u, u, u, u}, + {o, o, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + }, + }, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + targetLevel := Level(tt.tmsID) //nolint:gosec // G115 + segments, classification := tt.ix.lineTracePolygonEdges(tt.polygon, tt.tmsID, tt.buffer) + tt.ix.fillTiles(targetLevel, 0, 0, true, segments, classification, tt.polygon) + + got := classificationGrid(classification, targetLevel) + for l := Level(0); l <= targetLevel; l++ { + t.Logf("level %d: %v", l, got[l]) + } + assert.Equal(t, tt.want, got) + }) + } +} diff --git a/pointindex/pointindex.go b/pointindex/pointindex.go index 52f7f39..32121c9 100644 --- a/pointindex/pointindex.go +++ b/pointindex/pointindex.go @@ -11,6 +11,7 @@ import ( "github.com/go-spatial/geom/encoding/wkt" "golang.org/x/exp/maps" //nolint:exptostd + "github.com/pdok/texel/geomhelp" "github.com/pdok/texel/mathhelp" "github.com/pdok/texel/morton" "github.com/pdok/texel/tms20" @@ -42,6 +43,14 @@ type Quadrant struct { intCentroid intgeom.Point } +func (q *Quadrant) Extent() geom.Extent { + return q.intExtent.ToGeomExtent() +} + +func (q *Quadrant) Coords() (uint, uint) { + return morton.FromZ(q.z) +} + // PointIndex is a pointcloud annex quadtree to enable snapping lines to a grid accounting for those points. // Quadrants: // @@ -66,20 +75,25 @@ type PointIndex struct { Quadrant 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 + 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 + tilePixels uint + internalPixels uint } -type Level = uint -type Q = int // quadrant index (0, 1, 2 or 3) +type ( + Level = uint + Q = int // quadrant index (0, 1, 2 or 3) +) func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID) (*PointIndex, error) { // assuming IsQuadTree was tested before rootTM := tileMatrixSet.TileMatrices[0] - levelDiff := uint(math.Log2(float64(rootTM.TileWidth))) + uint(math.Log2(float64(VectorTileInternalPixelResolution))) + tilePixels := rootTM.TileWidth + levelDiff := uint(math.Log2(float64(tilePixels))) + uint(math.Log2(float64(VectorTileInternalPixelResolution))) //nolint:gosec // G115 deepestLevel := uint(deepestTMID) + levelDiff bottomLeft, topRight, err := tileMatrixSet.MatrixBoundingBox(0) @@ -95,8 +109,10 @@ func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID intExtent: intExtent, z: 0, }, - deepestLevel: deepestLevel, - deepestSize: deepestSize, + tilePixels: tilePixels, + internalPixels: VectorTileInternalPixelResolution, + deepestLevel: deepestLevel, + deepestSize: deepestSize, //nolint:gosec // G115 deepestRes: intExtent.XSpan() / int64(deepestSize), quadrants: make(map[Level]map[morton.Z]Quadrant, deepestLevel+1), @@ -108,6 +124,35 @@ func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID return &ix, nil } +func Factory(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID) func() *PointIndex { + factory := func() *PointIndex { + ix, err := FromTileMatrixSet(tileMatrixSet, deepestTMID) + if err != nil { + panic(err) + } + return ix + } + return factory +} + +// Insert all points of a geometry in the pointindex. Multi-geometries not supported +func (ix *PointIndex) InsertGeometry(geometry geom.Geometry) error { + pointSlice := geomhelp.PointSlice(geometry) + // Initialize the map + for level := range ix.deepestLevel + 1 { + if ix.quadrants[level] == nil { + ix.quadrants[level] = make(map[morton.Z]Quadrant, len(pointSlice)) + } + } + // Insert points + for _, point := range pointSlice { + if err := ix.InsertPoint(point); err != nil { + return err + } + } + return nil +} + // InsertPolygon inserts all points from a Polygon func (ix *PointIndex) InsertPolygon(polygon geom.Polygon) error { // initialize the quadrants map diff --git a/pointindex/pointindex_test.go b/pointindex/pointindex_test.go index 40927ea..6a219f9 100644 --- a/pointindex/pointindex_test.go +++ b/pointindex/pointindex_test.go @@ -148,8 +148,10 @@ 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}), }, - deepestLevel: 0, - deepestSize: mathhelp.Pow2(0), + deepestLevel: 0, + deepestSize: mathhelp.Pow2(0), + tilePixels: 256, + internalPixels: 16, //nolint:gosec // G115 deepestRes: intgeom.FromGeomOrd(1.0) / intgeom.M(mathhelp.Pow2(0)), quadrants: map[Level]map[morton.Z]Quadrant{0: {0: Quadrant{ @@ -167,8 +169,10 @@ 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}), }, - deepestLevel: 1, - deepestSize: mathhelp.Pow2(1), + deepestLevel: 1, + deepestSize: mathhelp.Pow2(1), + tilePixels: 256, + internalPixels: 16, //nolint:gosec // G115 deepestRes: intgeom.FromGeomOrd(1.0) / intgeom.M(mathhelp.Pow2(1)), quadrants: map[Level]map[morton.Z]Quadrant{ @@ -193,8 +197,10 @@ 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}), }, - deepestLevel: 3, - deepestSize: mathhelp.Pow2(3), + deepestLevel: 3, + deepestSize: mathhelp.Pow2(3), + tilePixels: 256, + internalPixels: 16, //nolint:gosec // G115 deepestRes: intgeom.FromGeomOrd(4.0) / intgeom.M(mathhelp.Pow2(3)), quadrants: map[Level]map[morton.Z]Quadrant{ @@ -230,8 +236,10 @@ 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}), }, - deepestLevel: 5, - deepestSize: mathhelp.Pow2(5), + deepestLevel: 5, + deepestSize: mathhelp.Pow2(5), + tilePixels: 256, + internalPixels: 16, //nolint:gosec // G115 deepestRes: intgeom.FromGeomOrd(16.0) / intgeom.M(mathhelp.Pow2(5)), quadrants: map[Level]map[morton.Z]Quadrant{ @@ -473,7 +481,7 @@ func TestPointIndex_SnapClosestPoints(t *testing.T) { t.Run(tt.name, func(t *testing.T) { ix := tt.ix poly := tt.poly - err := ix.InsertPolygon(poly) + err := ix.InsertGeometry(poly) require.NoError(t, err) levels := tt.levels if levels == nil { @@ -509,6 +517,17 @@ func TestPointIndex_lineIntersects(t *testing.T) { }, want: false, }, + { + name: "diagonally opposed line between non-inclusive points", + // This test is fragile and depends on floating-point rounding + extent: intgeom.Extent{ + 0o0000000, 0o0000000, 10000000, 10000000, + }, + line: intgeom.Line{ + {0o0000000, 10000000}, {10000000, 0o0000000}, + }, + want: false, // TODO This is an undesired outcome + }, } for _, tt := range tests { t.Run(tt.name, func(t *testing.T) { @@ -532,15 +551,37 @@ func newSimplePointIndex(deepestLevel Level, cellSize float64) *PointIndex { deepestLevel: deepestLevel, deepestSize: deepestSize, //nolint:gosec // G115 - 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), + 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), + tilePixels: 256, + internalPixels: 16, } _, ix.intCentroid = ix.getQuadrantExtentAndCentroid(0, 0, 0, ix.intExtent) return &ix } +func newPointIndex(deepestTMID tms20.TMID, tilePixels, internalPixels uint, cellSize float64) *PointIndex { + ix := PointIndex{ + hitOnce: make(map[morton.Z]map[intgeom.Point][]int, 0), + hitMultiple: make(map[morton.Z]map[intgeom.Point][]int, 0), + tilePixels: tilePixels, + internalPixels: internalPixels, + } + deepestLevel := ix.InternalPixelLevelFromTmsID(deepestTMID) + ix.deepestLevel = deepestLevel + ix.deepestSize = mathhelp.Pow2(deepestLevel) + span := cellSize * float64(ix.deepestSize) + intExtent := intgeom.Extent{0.0, 0.0, intgeom.FromGeomOrd(span), intgeom.FromGeomOrd(span)} + ix.Quadrant = Quadrant{ + intExtent: intExtent, + } + ix.deepestRes = intExtent.XSpan() / int64(ix.deepestSize) //nolint:gosec // G115 + ix.quadrants = make(map[Level]map[morton.Z]Quadrant, deepestLevel+1) + return &ix +} + func loadEmbeddedTileMatrixSet(t *testing.T, tmsID string) tms20.TileMatrixSet { t.Helper() diff --git a/processing/gpkg/gpkg.go b/processing/gpkg/gpkg.go index 8c4328f..b6af9c0 100644 --- a/processing/gpkg/gpkg.go +++ b/processing/gpkg/gpkg.go @@ -1,29 +1,16 @@ package gpkg +// This file defines general geopackage-related infrastructure shared by both +// the `snap` and `mvt` flows. + import ( "fmt" "log" "strings" - "time" - "github.com/go-spatial/geom" "github.com/go-spatial/geom/encoding/gpkg" - "github.com/pdok/texel/processing" ) -type featureGPKG struct { - columns []any - geometry geom.Geometry -} - -func (f featureGPKG) Columns() []any { - return f.columns -} - -func (f featureGPKG) Geometry() geom.Geometry { - return f.geometry -} - type column struct { cid int name string @@ -33,6 +20,7 @@ type column struct { pk int } +// Abstraction of geometry table in a gpkg type Table struct { Name string columns []column @@ -41,7 +29,11 @@ type Table struct { srs gpkg.SpatialReferenceSystem } -// geometryTypeFromString returns the numeric value of a gometry string +// Name of associated table of a geometry table. +func (t Table) EncodedName() string { + return t.Name + "_encoded" +} + func geometryTypeFromString(geometrytype string) gpkg.GeometryType { switch strings.ToUpper(geometrytype) { case "GEOMETRY": @@ -65,90 +57,21 @@ func geometryTypeFromString(geometrytype string) gpkg.GeometryType { } } -//nolint:recvcheck // TODO double check this -type SourceGeopackage struct { - Table Table - handle *gpkg.Handle -} +// Custom wrapper for gpkg.Handle. To be embedded in other structs +type geopackageHandle struct{ handle *gpkg.Handle } -func (source *SourceGeopackage) Init(file string) { - source.handle = openGeopackage(file) +func (gH *geopackageHandle) Init(file string) { + gH.handle = openGeopackage(file) } -func (source SourceGeopackage) Close() { - source.handle.Close() +func (gH *geopackageHandle) Close() { + gH.handle.Close() } -//nolint:funlen,cyclop -func (source SourceGeopackage) ReadFeatures(features chan<- processing.Feature) { - - rows, err := source.handle.Query(source.Table.selectSQL()) - if err != nil { - log.Fatalf("err during closing rows: %s", err) - } - - cols, err := rows.Columns() - if err != nil { - log.Fatalf("error reading the columns: %s", err) - } - - for rows.Next() { - vals := make([]any, len(cols)) - valPtrs := make([]any, len(cols)) - for i := range cols { - valPtrs[i] = &vals[i] - } - - if err = rows.Scan(valPtrs...); err != nil { - log.Fatalf("err reading row values: %v", err) - } - var f featureGPKG - var c []any - - for i, colName := range cols { - switch colName { - case source.Table.gcolumn: - wkbgeom, err := gpkg.DecodeGeometry(vals[i].([]byte)) - if err != nil { - log.Fatalf("error decoding the geometry: %s", err) - } - f.geometry = wkbgeom.Geometry - default: - switch v := vals[i].(type) { - case []uint8: - asBytes := make([]byte, len(v)) - copy(asBytes, v) - c = append(c, string(asBytes)) - case int64: - c = append(c, v) - case float64: - c = append(c, v) - case time.Time: - c = append(c, v.Format(time.RFC3339)) - case string: - c = append(c, v) - case nil: - c = append(c, v) - default: - log.Fatalf("unexpected type for sqlite column data: %v: %T", cols[i], v) - } - } - f.columns = c - } - ff := &f - features <- ff - } - err = rows.Err() - if err != nil { - log.Fatal(err) - } - close(features) - defer rows.Close() -} - -func (source SourceGeopackage) GetTableInfo() []Table { +// Generate []Table slice for looping over. +func (gH *geopackageHandle) GetTableInfo() []Table { query := `SELECT table_name, column_name, geometry_type_name, srs_id FROM gpkg_geometry_columns;` - rows, err := source.handle.Query(query) + rows, err := gH.handle.Query(query) if err != nil { log.Fatalf("error during closing rows: %v - %v", query, err) } @@ -163,9 +86,9 @@ func (source SourceGeopackage) GetTableInfo() []Table { log.Fatalf("error retrieving the source table information: %s", err) } - t.columns = getTableColumns(source.handle, t.Name) + t.columns = gH.getTableColumns(t.Name) t.gtype = geometryTypeFromString(gtype) - t.srs = getSpatialReferenceSystem(source.handle, srsID) + t.srs = gH.getSpatialReferenceSystem(srsID) tables = append(tables, t) } @@ -173,106 +96,6 @@ func (source SourceGeopackage) GetTableInfo() []Table { return tables } -type TargetGeopackage struct { - Table Table - pagesize int - handle *gpkg.Handle -} - -func (target *TargetGeopackage) Init(file string, pagesize int) { - target.pagesize = pagesize - target.handle = openGeopackage(file) -} - -func (target *TargetGeopackage) Close() { - target.handle.Close() -} - -func (target *TargetGeopackage) CreateTables(tables []Table) error { - for _, table := range tables { - err := target.handle.UpdateSRS(table.srs) - if err != nil { - return err - } - - err = buildTable(target.handle, table) - if err != nil { - return err - } - } - return nil -} - -func (target *TargetGeopackage) WriteFeatures(inFeatures <-chan processing.Feature) { - var features []processing.Feature - - for { - feature, hasMore := <-inFeatures - if !hasMore { - target.writeFeatures(features) - break - } - features = append(features, feature) - - if len(features)%target.pagesize == 0 { - target.writeFeatures(features) - features = nil - } - } -} - -func (target *TargetGeopackage) writeFeatures(features []processing.Feature) { - tx, err := target.handle.Begin() - if err != nil { - log.Fatalf("Could not start a transaction: %s", err) - } - - stmt, err := tx.Prepare(target.Table.insertSQL()) - if err != nil { - log.Fatalf("Could not prepare a statement: %s", err) - } - - var ext *geom.Extent - - for _, f := range features { - //nolint:gosec // G115 - sb, err := gpkg.NewBinary(int32(target.Table.srs.ID), f.Geometry()) - if err != nil { - log.Fatalf("Could not create a binary geometry: %s", err) - } - - data := f.Columns() - data = append(data, sb) - - _, err = stmt.Exec(data...) - if err != nil { - var fid any = "unknown" - if len(data) > 0 { - fid = data[0] - } - log.Fatalf("Could not get a result summary from the prepared statement for fid %s: %s", fid, err) - } - - if ext == nil { - ext, err = geom.NewExtentFromGeometry(f.Geometry()) - if err != nil { - ext = nil - log.Println("Failed to create new extent:", err) - continue - } - } else { - _ = ext.AddGeometry(f.Geometry()) - } - } - stmt.Close() - _ = tx.Commit() - - err = target.handle.UpdateGeometryExtent(target.Table.Name, ext) - if err != nil { - log.Fatalln("Failed to update new extent:", err) - } -} - func openGeopackage(file string) *gpkg.Handle { handle, err := gpkg.Open(file) if err != nil { @@ -281,60 +104,12 @@ func openGeopackage(file string) *gpkg.Handle { return handle } -// createSQL creates a CREATE statement on the given table and column information -// used for creating feature tables in the target Geopackage -func (t Table) createSQL() string { - create := fmt.Sprintf(`CREATE TABLE IF NOT EXISTS "%v"`, t.Name) - var columnparts []string //nolint:prealloc - for _, column := range t.columns { - columnpart := column.name + ` ` + column.ctype - if column.notnull == 1 { - columnpart += ` NOT NULL` - } - if column.pk == 1 { - columnpart += ` PRIMARY KEY` - } - - columnparts = append(columnparts, columnpart) - } - - query := create + `(` + strings.Join(columnparts, `, `) + `);` - return query -} - -// selectSQL build a SELECT statement based on the table and columns -// used for reading the source features -func (t Table) selectSQL() string { - var csql []string //nolint:prealloc - for _, c := range t.columns { - csql = append(csql, c.name) - } - query := `SELECT ` + strings.Join(csql, `,`) + ` FROM "` + t.Name + `";` - return query -} - -// insertSQL used for writing the features -// build the INSERT statement based on the table and columns -func (t Table) insertSQL() string { - var csql, vsql []string - for _, c := range t.columns { - if c.name != t.gcolumn { - csql = append(csql, c.name) - vsql = append(vsql, `?`) - } - } - csql = append(csql, t.gcolumn) - vsql = append(vsql, `?`) - query := `INSERT INTO "` + t.Name + `"(` + strings.Join(csql, `,`) + `) VALUES(` + strings.Join(vsql, `,`) + `)` - return query -} - // getSpatialReferenceSystem extracts this based on the given SRS id -func getSpatialReferenceSystem(h *gpkg.Handle, id int) gpkg.SpatialReferenceSystem { +func (gH *geopackageHandle) getSpatialReferenceSystem(id int) gpkg.SpatialReferenceSystem { var srs gpkg.SpatialReferenceSystem query := `SELECT srs_name, srs_id, organization, organization_coordsys_id, definition, description FROM gpkg_spatial_ref_sys WHERE srs_id = %v;` - row := h.QueryRow(fmt.Sprintf(query, id)) + row := gH.handle.QueryRow(fmt.Sprintf(query, id)) var description *string _ = row.Scan(&srs.Name, &srs.ID, &srs.Organization, &srs.OrganizationCoordsysID, &srs.Definition, &description) if description != nil { @@ -345,11 +120,10 @@ func getSpatialReferenceSystem(h *gpkg.Handle, id int) gpkg.SpatialReferenceSyst } // getTableColumns collects the column information of a given table -func getTableColumns(h *gpkg.Handle, table string) []column { +func (gH *geopackageHandle) getTableColumns(table string) []column { var columns []column query := `PRAGMA table_info('%v');` - rows, err := h.Query(fmt.Sprintf(query, table)) - + rows, err := gH.handle.Query(fmt.Sprintf(query, table)) if err != nil { log.Fatalf("err during closing rows: %v - %v", query, err) } @@ -365,30 +139,3 @@ func getTableColumns(h *gpkg.Handle, table string) []column { defer rows.Close() return columns } - -// buildTable creates a given destination table with the necessary gpkg_ information -func buildTable(h *gpkg.Handle, t Table) error { - query := t.createSQL() - _, err := h.Exec(query) - if err != nil { - log.Fatalf("error building table in target GeoPackage: %s", err) - } - - err = h.AddGeometryTable(gpkg.TableDescription{ - Name: t.Name, - ShortName: t.Name, - Description: t.Name, - GeometryField: t.gcolumn, - GeometryType: t.gtype, - //nolint:gosec // G115 - SRS: int32(t.srs.ID), - // - Z: gpkg.Prohibited, - M: gpkg.Prohibited, - }) - if err != nil { - log.Println("error adding geometry table in target GeoPackage:", err) - return err - } - return nil -} diff --git a/processing/gpkg/gpkgmvt.go b/processing/gpkg/gpkgmvt.go new file mode 100644 index 0000000..741e4ea --- /dev/null +++ b/processing/gpkg/gpkgmvt.go @@ -0,0 +1,235 @@ +package gpkg + +// Specializes the gpkg.go functionaltity for the `texel mvt` command + +import ( + "encoding/binary" + "fmt" + "strings" + "time" + + "github.com/pdok/texel/processing" + "github.com/pdok/texel/tile" +) + +// Batch value for reading features from the geopackage table +const maxSQLiteBatchParams = 500 + +type MVTSourceGeopackage struct { + Table Table + geopackageHandle +} + +// GetAttributesForFeatures fetches the attributes (i.e. all non-geometry, non-fid columns) +// for the given feature IDs, batching the query to respect SQLite's parameter limits. +// The fid is assumed to be the table's first column. +func (source MVTSourceGeopackage) GetAttributesForFeatures(featureIDs []int64) (tile.InternalAttributeTable, error) { + attributes := make(tile.InternalAttributeTable, len(featureIDs)) + fidColumn := source.Table.columns[0].name + + for start := 0; start < len(featureIDs); start += maxSQLiteBatchParams { + batch := featureIDs[start:min(start+maxSQLiteBatchParams, len(featureIDs))] + + batchAsAny := make([]any, len(batch)) + for i, fid := range batch { + batchAsAny[i] = fid + } + + if err := source.queryAttributesBatch(fidColumn, batchAsAny, attributes); err != nil { + return nil, err + } + } + + return attributes, nil +} + +// selectByFIDsSQL builds a SELECT of all the table's columns for a batch of feature IDs. +// The fid is assumed to be the table's first column. +func (t Table) selectByFIDsSQL(batchSize int) string { + var csql []string //nolint:prealloc + for _, c := range t.columns { + csql = append(csql, c.name) + } + placeholders := strings.TrimSuffix(strings.Repeat("?,", batchSize), ",") + fidColumn := t.columns[0].name + return `SELECT ` + strings.Join(csql, `,`) + ` FROM "` + t.Name + `" WHERE ` + fidColumn + ` IN (` + placeholders + `)` +} + +func cleanData(v any) (any, bool) { + switch v := v.(type) { + case []uint8: + asBytes := make([]byte, len(v)) + copy(asBytes, v) + return v, true + case int64: + return v, true + case float64: + return v, true + case time.Time: + return v.Format(time.RFC3339), true + case string: + return v, true + case nil: + return nil, true + default: + return nil, false + } +} + +// attributeColumnNames lists all columns except fid and geometry column. +// Assumptions: fid column has index 0 +func (t Table) attributeColumnNames() []string { + names := make([]string, 0, len(t.columns)) + for i, c := range t.columns { + if i == 0 || c.name == t.gcolumn { + continue // skip the fid and geometry columns, they are not attributes + } + names = append(names, c.name) + } + return names +} + +// AttributeColumnNames satisfies processing.MVTSource by delegating to the +// current Table. +func (source MVTSourceGeopackage) AttributeColumnNames() []string { + return source.Table.attributeColumnNames() +} + +// selectEncodedSQL builds the SELECT for fetching the encoded features of a single tile +func (t Table) selectEncodedSQL() string { + return `SELECT feature_id, geometry_type, data FROM "` + t.EncodedName() + `" WHERE tile_x = ? AND tile_y = ? AND tile_z = ?` +} + +// selectDistinctTilesSQL builds the SELECT for enumerating every tile that has at +// least one encoded feature. +func (t Table) selectDistinctTilesSQL() string { + return `SELECT DISTINCT tile_x, tile_y FROM "` + t.EncodedName() + `" WHERE tile_z = ? ORDER BY tile_x, tile_y` +} + +// GetFeaturesForTile returns every encoded feature stored for the given tile, +// i.e. the equivalent of the tile-features table's rows for (tileX, tileY). +func (source MVTSourceGeopackage) GetFeaturesForTile(tileX, tileY, tileZ uint) ([]tile.EncodedFeatureRow, error) { + rows, err := source.handle.Query(source.Table.selectEncodedSQL(), tileX, tileY, tileZ) + if err != nil { + return nil, fmt.Errorf("querying encoded features for tile (%d,%d,%d): %w", tileX, tileY, tileZ, err) + } + defer rows.Close() + + var features []tile.EncodedFeatureRow + for rows.Next() { + var featureID int64 + var geometryType int32 + var data []byte + if err := rows.Scan(&featureID, &geometryType, &data); err != nil { + return nil, fmt.Errorf("scanning encoded feature row for tile (%d,%d): %w", tileX, tileY, err) + } + features = append(features, tile.EncodedFeatureRow{ + FeatureID: featureID, + Geom: tile.EncodedGeometry{ + Encoding: deserializeFromBytes(data), + GeometryType: geometryType, + XTile: tileX, + YTile: tileY, + }, + }) + } + if err := rows.Err(); err != nil { + return nil, fmt.Errorf("iterating encoded feature rows for tile (%d,%d): %w", tileX, tileY, err) + } + return features, nil +} + +// deserializeFromBytes is the inverse of serializeToBytes: it turns a BLOB of +// little-endian uint32s back into the encoded geometry command/coordinate stream. +func deserializeFromBytes(data []byte) []uint32 { + intslice := make([]uint32, len(data)/4) + for i := range intslice { + intslice[i] = binary.LittleEndian.Uint32(data[i*4:]) + } + return intslice +} + +// ListTiles returns every distinct tile present in the "