Second attempt at ST_AsMVT and ST_AsMVTGeom

docs/functions.md

@@ -12,6 +12,7 @@
 | [`ST_Area_Spheroid`](#st_area_spheroid) | Returns the area of a geometry in meters, using an ellipsoidal model of the earth |
 | [`ST_AsGeoJSON`](#st_asgeojson) | Returns the geometry as a GeoJSON fragment |
 | [`ST_AsHEXWKB`](#st_ashexwkb) | Returns the geometry as a HEXWKB string |
+| [`ST_AsMVTGeom`](#st_asmvtgeom) | Transform and clip geometry to a tile boundary |
 | [`ST_AsSVG`](#st_assvg) | Convert the geometry into a SVG fragment or path |
 | [`ST_AsText`](#st_astext) | Returns the geometry as a WKT string |
 | [`ST_AsWKB`](#st_aswkb) | Returns the geometry as a WKB (Well-Known-Binary) blob |
@@ -88,6 +89,7 @@
 | [`ST_MakeBox2D`](#st_makebox2d) | Create a BOX2D from two POINT geometries |
 | [`ST_MakeEnvelope`](#st_makeenvelope) | Create a rectangular polygon from min/max coordinates |
 | [`ST_MakeLine`](#st_makeline) | Create a LINESTRING from a list of POINT geometries |
+| [`ST_MakePoint`](#st_makepoint) | Creates a GEOMETRY point from an pair of floating point numbers. |
 | [`ST_MakePolygon`](#st_makepolygon) | Create a POLYGON from a LINESTRING shell |
 | [`ST_MakeValid`](#st_makevalid) | Returns a valid representation of the geometry |
 | [`ST_MaximumInscribedCircle`](#st_maximuminscribedcircle) | Returns the maximum inscribed circle of the input geometry, optionally with a tolerance. |
@@ -144,6 +146,7 @@
 
 | Function | Summary |
 | --- | --- |
+| [`ST_AsMVT`](#st_asmvt) | Make a Mapbox Vector Tile from a set of geometries and properties |
 | [`ST_CoverageInvalidEdges_Agg`](#st_coverageinvalidedges_agg) | Returns the invalid edges of a coverage geometry |
 | [`ST_CoverageSimplify_Agg`](#st_coveragesimplify_agg) | Simplifies a set of geometries while maintaining coverage |
 | [`ST_CoverageUnion_Agg`](#st_coverageunion_agg) | Unions a set of geometries while maintaining coverage |
@@ -402,6 +405,26 @@ SELECT ST_AsHexWKB('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))'::geometry);
 
 ----
 
+### ST_AsMVTGeom
+
+
+#### Signatures
+
+```sql
+GEOMETRY ST_AsMVTGeom (geom GEOMETRY, bounds BOX_2D, extent BIGINT, buffer BIGINT, clip_geom BOOLEAN)
+GEOMETRY ST_AsMVTGeom (geom GEOMETRY, bounds BOX_2D, extent BIGINT, buffer BIGINT)
+GEOMETRY ST_AsMVTGeom (geom GEOMETRY, bounds BOX_2D, extent BIGINT)
+GEOMETRY ST_AsMVTGeom (geom GEOMETRY, bounds BOX_2D)
+```
+
+#### Description
+
+Transform and clip geometry to a tile boundary
+
+See "ST_AsMVT" for more details
+
+----
+
 ### ST_AsSVG
 
 
@@ -1937,6 +1960,35 @@ LINESTRING(0 0, 1 1)
 
 ----
 
+### ST_MakePoint
+
+
+#### Signatures
+
+```sql
+POINT_2D ST_MakePoint (x DOUBLE, y DOUBLE)
+POINT_3D ST_MakePoint (x DOUBLE, y DOUBLE, z DOUBLE)
+POINT_4D ST_MakePoint (x DOUBLE, y DOUBLE, z DOUBLE, m DOUBLE)
+```
+
+#### Description
+
+Creates a GEOMETRY point from an pair of floating point numbers.
+
+For geodetic coordinate systems, x is typically the longitude value and y is the latitude value.
+
+Note that ST_Point is equivalent. ST_MakePoint is provided for PostGIS compatibility.
+
+#### Example
+
+```sql
+SELECT ST_AsText(ST_MakePoint(143.3, -24.2));
+----
+POINT (143.3 -24.2)
+```
+
+----
+
 ### ST_MakePolygon
 
 
@@ -3015,6 +3067,66 @@ SELECT ST_ZMin(ST_Point(1, 2, 3))
 
 ## Aggregate Functions
 
+### ST_AsMVT
+
+
+#### Signatures
+
+```sql
+BLOB ST_AsMVT (col0 ANY)
+BLOB ST_AsMVT (col0 ANY, col1 VARCHAR)
+BLOB ST_AsMVT (col0 ANY, col1 VARCHAR, col2 INTEGER)
+BLOB ST_AsMVT (col0 ANY, col1 VARCHAR, col2 INTEGER, col3 VARCHAR)
+BLOB ST_AsMVT (col0 ANY, col1 VARCHAR, col2 INTEGER, col3 VARCHAR, col4 VARCHAR)
+```
+
+#### Description
+
+Make a Mapbox Vector Tile from a set of geometries and properties
+The function takes as input a row type (STRUCT) containing a geometry column and any number of property columns.
+It returns a single binary BLOB containing the Mapbox Vector Tile.
+
+The function has the following signature:
+
+`ST_AsMVT(row STRUCT, layer_name VARCHAR DEFAULT 'layer', extent INTEGER DEFAULT 4096, geom_column_name VARCHAR DEFAULT NULL, feature_id_column_name VARCHAR DEFAULT NULL) -> BLOB`
+
+- The first argument is a struct containing the geometry and properties.
+- The second argument is the name of the layer in the vector tile. This argument is optional and defaults to 'layer'.
+- The third argument is the extent of the tile. This argument is optional and defaults to 4096.
+- The fourth argument is the name of the geometry column in the input row. This argument is optional. If not provided, the first geometry column in the input row will be used. If multiple geometry columns are present, an error will be raised.
+- The fifth argument is the name of the feature id column in the input row. This argument is optional. If provided, the values in this column will be used as feature ids in the vector tile. The column must be of type INTEGER or BIGINT. If set to negative or NULL, a feature id will not be assigned to the corresponding feature.
+
+The input struct must contain exactly one geometry column of type GEOMETRY. It can contain any number of property columns of types VARCHAR, FLOAT, DOUBLE, INTEGER, BIGINT, or BOOLEAN.
+
+Example:
+```sql
+SELECT ST_AsMVT({'geom': geom, 'id': id, 'name': name}, 'cities', 4096, 'geom', 'id') AS tile
+FROM cities;
+ ```
+
+This example creates a vector tile named 'cities' with an extent of 4096 from the 'cities' table, using 'geom' as the geometry column and 'id' as the feature id column.
+
+However, you probably want to use the ST_AsMVTGeom function to first transform and clip your geometries to the tile extent.
+The following example assumes the geometry is in WebMercator ("EPSG:3857") coordinates.
+Replace `{z}`, `{x}`, and `{y}` with the appropriate tile coordinates, `{your table}` with your table name, and `{tile_path}` with the path to write the tile to.
+
+```sql
+COPY (
+    SELECT ST_AsMVT({{
+        "geometry": ST_AsMVTGeom(
+            geometry,
+            ST_Extent(ST_TileEnvelope({z}, {x}, {y})),
+            4096,
+            256,
+            false
+        )
+    }})
+    FROM {your table} WHERE ST_Intersects(geometry, ST_TileEnvelope({z}, {x}, {y}))
+) to {tile_path} (FORMAT 'BLOB');
+```
+
+----
+
 ### ST_CoverageInvalidEdges_Agg
 
 

src/sgl/sgl.cpp

@@ -2976,7 +2976,7 @@ point_in_polygon_result prepared_geometry::contains(const vertex_xy &vert) const
 
 			const auto end = math::min(node_end, levl_end);
 
-			if(stack[depth] != end) {
+			if (stack[depth] != end) {
 				// Go sideways!
 				stack[depth]++;
 				break;
@@ -3367,7 +3367,7 @@ static bool try_get_prepared_distance_lines(const prepared_geometry &lhs, const
 
 	if (found_any) {
 		distance = std::sqrt(min_dist); // Convert squared distance to actual distance
-		return true; // We found a distance
+		return true;                    // We found a distance
 	}
 	return false; // No distance found
 }
@@ -3381,7 +3381,6 @@ bool prepared_geometry::try_get_distance(const prepared_geometry &other, double
 // WKT Parsing
 //======================================================================================================================
 
-
 namespace sgl {
 
 namespace {

src/spatial/modules/geos/geos_geometry.hpp

@@ -51,6 +51,13 @@ class GeosGeometry {
 	GeosGeometry get_voronoi_diagram() const;
 	GeosGeometry get_built_area() const;
 	GeosGeometry get_noded() const;
+	GeosGeometry get_clipped(double xmin, double ymin, double xmax, double ymax) const;
+
+	// matrix format: [a, b, c, d, e, f]
+	//	x' = a*x + b*y + e
+	//	y' = c*x + d*y + f
+	GeosGeometry get_transformed(const double matrix[6]) const;
+	GeosGeometry get_gridded(double grid_size) const;
 
 	bool contains(const GeosGeometry &other) const;
 	bool covers(const GeosGeometry &other) const;
@@ -67,6 +74,7 @@ class GeosGeometry {
 	double distance_to(const GeosGeometry &other) const;
 
 	void normalize_in_place() const;
+	void orient_polygons(bool ext_cw);
 
 	GeosGeometry get_difference(const GeosGeometry &other) const;
 	GeosGeometry get_intersection(const GeosGeometry &other) const;
@@ -94,6 +102,10 @@ class GeosGeometry {
 
 	PreparedGeosGeometry get_prepared() const;
 
+	void get_extent(double &xmin, double &ymin, double &xmax, double &ymax) const {
+		GEOSGeom_getExtent_r(handle, geom, &xmin, &ymin, &xmax, &ymax);
+	}
+
 private:
 	GEOSContextHandle_t handle;
 	GEOSGeometry *geom;
@@ -332,6 +344,34 @@ inline GeosGeometry GeosGeometry::get_noded() const {
 	return GeosGeometry(handle, GEOSNode_r(handle, geom));
 }
 
+inline GeosGeometry GeosGeometry::get_clipped(double xmin, double ymin, double xmax, double ymax) const {
+	return GeosGeometry(handle, GEOSClipByRect_r(handle, geom, xmin, ymin, xmax, ymax));
+}
+
+inline GeosGeometry GeosGeometry::get_transformed(const double matrix[6]) const {
+	// x' = a*x + b*y + e
+	// y' = c*x + d*y + f
+	return GeosGeometry(handle, GEOSGeom_transformXY_r(
+	                                handle, geom,
+	                                [](double *x_ptr, double *y_ptr, void *data) -> int {
+		                                const auto m = static_cast<const double *>(data);
+		                                const auto &x = *x_ptr;
+		                                const auto &y = *y_ptr;
+
+		                                const auto new_x = m[0] * x + m[1] * y + m[4];
+		                                const auto new_y = m[2] * x + m[3] * y + m[5];
+
+		                                *x_ptr = new_x;
+		                                *y_ptr = new_y;
+		                                return 1;
+	                                },
+	                                const_cast<double *>(matrix)));
+}
+
+inline GeosGeometry GeosGeometry::get_gridded(double grid_size) const {
+	return GeosGeometry(handle, GEOSGeom_setPrecision_r(handle, geom, grid_size, GEOS_PREC_NO_TOPO));
+}
+
 inline GeosGeometry GeosGeometry::get_maximum_inscribed_circle() const {
 	double xmin = 0;
 	double ymin = 0;
@@ -411,6 +451,10 @@ inline void GeosGeometry::normalize_in_place() const {
 	GEOSNormalize_r(handle, geom);
 }
 
+inline void GeosGeometry::orient_polygons(bool ext_cw) {
+	GEOSOrientPolygons_r(handle, geom, ext_cw ? 1 : 0);
+}
+
 inline GeosGeometry GeosGeometry::get_difference(const GeosGeometry &other) const {
 	return GeosGeometry(handle, GEOSDifference_r(handle, geom, other.geom));
 }