From 7d1274e9520a9b98d340ba9bb814153a3f8c4ed5 Mon Sep 17 00:00:00 2001 From: tidwall Date: Wed, 1 Jan 2025 04:55:40 -0700 Subject: [PATCH] Add support for GeoBIN See docs/API.md --- README.md | 7 +- docs/API.md | 166 +++++++++++++++- docs/GEOBIN.md | 137 +++++++++++++ tests/test_geobin.c | 173 +++++++++++++++++ tests/test_rect.c | 63 ++++++ tests/tests.h | 28 +++ tg.c | 464 ++++++++++++++++++++++++++++++++++++++++++-- tg.h | 6 + 8 files changed, 1020 insertions(+), 24 deletions(-) create mode 100644 docs/GEOBIN.md create mode 100644 tests/test_geobin.c diff --git a/README.md b/README.md index a49b324..2e062d2 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ I designed it for programs that need real-time geospatial, such as geofencing, m - Implements OGC [Simple Features](https://en.wikipedia.org/wiki/Simple_Features) including Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon, GeometryCollection. - Optimized [polygon indexing](docs/POLYGON_INDEXING.md) that introduces two new structures. -- Reads and writes [WKT](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry), [WKB](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry), and [GeoJSON](https://en.wikipedia.org/wiki/GeoJSON). +- Reads and writes [GeoJSON](https://en.wikipedia.org/wiki/GeoJSON), [WKT](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry), [WKB](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry), and [GeoBIN](docs/GEOBIN.md). - Provides a purely functional [API](docs/API.md) that is reentrant and thread-safe. - Spatial predicates including "intersects", "covers", "touches", "equals", etc. - Compiles to Webassembly using Emscripten @@ -30,8 +30,7 @@ It's a non-goal for TG to be a full GIS library. Consider [GEOS](https://libgeos TG uses [entirely new](docs/POLYGON_INDEXING.md) indexing structures that speed up [geometry predicates](docs/API.md#geometry-predicates). It can index more than 10GB per second of point data on modern hardware, while using less than 7% of additional memory, and can perform over 10 million point-in-polygon operations per second, even when using large polygons with over 10K points. -The following benchmark provides an example of the point-in-polygon performance -of TG when using a large polygon. In this case of Brazil, which has 39K points. +The following benchmark provides an example of the parsing, indexing, and point-in-polygon speed of TG vs GEOS when using a large polygon. In this case Brazil, which has 39K points. 
 Brazil              ops/sec    ns/op  points  hits       built      bytes
@@ -42,7 +41,7 @@ geos/none            29,708    33661   39914  3257   135.18 µs    958,104
 geos/prepared     7,885,512      127   39914  3257  2059.94 µs  3,055,496
-- "built": Column showing how much time the polygon and index took to construct. +- "built": Column showing how much time it took to construct the polygon and index. - "bytes": Column showing the final in-memory size of the polygon and index. - "none": No indexing was used. - "natural": Using TG [Natural](docs/POLYGON_INDEXING.md#natural) indexing diff --git a/docs/API.md b/docs/API.md index cbaf8d4..7e90096 100644 --- a/docs/API.md +++ b/docs/API.md @@ -182,6 +182,7 @@ Functions for accessing various information about geometries, such as getting th - [tg_geom_num_extra_coords()](#group___geometry_accessors_1ga5f85cf4c143703ee227e3c35accbde8b) - [tg_geom_memsize()](#group___geometry_accessors_1ga4931914f5170cce2949cfdd79e34ef63) - [tg_geom_search()](#group___geometry_accessors_1gad18cdb2a4ab1fa711dce821a0868ecd2) +- [tg_geom_fullrect()](#group___geometry_accessors_1gac1a077f09e247c022f09e48392b80051) @@ -224,9 +225,12 @@ Functions for parsing geometries from external data representations. It's recomm - [tg_parse_hexn()](#group___geometry_parsing_1ga2beb47ee201ddbd1a9a7a43c938c4396) - [tg_parse_hex_ix()](#group___geometry_parsing_1ga8718e134723418426e5df2b13acdf0c1) - [tg_parse_hexn_ix()](#group___geometry_parsing_1ga3d1f9eb85ff97812fed7cf7bf9e14bd4) +- [tg_parse_geobin()](#group___geometry_parsing_1ga4792bf4f319356fff3fa539cbd44b196) +- [tg_parse_geobin_ix()](#group___geometry_parsing_1gac0450996bcd71cdde81af451dd0e9571) - [tg_parse()](#group___geometry_parsing_1gaa9e5850bb2cc4eb442c227cd5ac738a1) - [tg_parse_ix()](#group___geometry_parsing_1gae0bfc62deb68979a46ed62facfee1280) - [tg_geom_error()](#group___geometry_parsing_1gae77b27ad34c2a215cc281647ab6dbc7e) +- [tg_geobin_fullrect()](#group___geometry_parsing_1gac77e3d8d51a7e66381627cb25853d80f) @@ -240,6 +244,7 @@ Functions for writing geometries as external data representations. - [tg_geom_wkt()](#group___geometry_writing_1ga047599bbe51886a6c3fbe35a6372d5d6) - [tg_geom_wkb()](#group___geometry_writing_1gac0938697f9270d1924c6746af68dd693) - [tg_geom_hex()](#group___geometry_writing_1ga269f61715302d48f144ffaf46e8c397b) +- [tg_geom_geobin()](#group___geometry_writing_1gab0f92319db4b8c61c62e0bd735d907b7) @@ -1629,6 +1634,32 @@ Iterates over all child geometries in geom that intersect rect + +## tg_geom_fullrect() +```c +int tg_geom_fullrect(const struct tg_geom *geom, double min[4], double max[4]); +``` +Returns the minimum bounding rectangle of a geometry on all dimensions. + +**Parameters** + +- **geom**: Input geometry +- **min**: min values, must have room for 4 dimensions +- **max**: max values, must have room for 4 dimensions + + + +**Return** + +- number of dimensions, or zero if invalid geom. + + +**See also** + +- [tg_geom_rect()](#group___geometry_accessors_1ga68d67f900b847ae08e6515a620f4f657) + + + ## tg_geom_equals() ```c @@ -2074,7 +2105,7 @@ Parse Well-known binary (WKB) using provided indexing option. ```c struct tg_geom *tg_parse_hex(const char *hex); ``` -Parse hex encoded Well-known binary (WKB). +Parse hex encoded Well-known binary (WKB) or GeoBIN. **Parameters** @@ -2103,7 +2134,7 @@ Parse hex encoded Well-known binary (WKB). ```c struct tg_geom *tg_parse_hexn(const char *hex, size_t len); ``` -Parse hex encoded Well-known binary (WKB) with an included data length. +Parse hex encoded Well-known binary (WKB) or GeoBIN with an included data length. **Parameters** @@ -2129,7 +2160,7 @@ Parse hex encoded Well-known binary (WKB) with an included data length. ```c struct tg_geom *tg_parse_hex_ix(const char *hex, enum tg_index ix); ``` -Parse hex encoded Well-known binary (WKB) using provided indexing option. +Parse hex encoded Well-known binary (WKB) or GeoBIN using provided indexing option. **Parameters** @@ -2156,7 +2187,7 @@ Parse hex encoded Well-known binary (WKB) using provided indexing option. ```c struct tg_geom *tg_parse_hexn_ix(const char *hex, size_t len, enum tg_index ix); ``` -Parse hex encoded Well-known binary (WKB) using provided indexing option. +Parse hex encoded Well-known binary (WKB) or GeoBIN using provided indexing option. **Parameters** @@ -2179,6 +2210,62 @@ Parse hex encoded Well-known binary (WKB) using provided indexing option. + +## tg_parse_geobin() +```c +struct tg_geom *tg_parse_geobin(const uint8_t *geobin, size_t len); +``` +Parse GeoBIN binary using provided indexing option. + +**Parameters** + +- **geobin**: GeoBIN data +- **len**: Length of data +- **ix**: Indexing option, e.g. TG_NONE, TG_NATURAL, TG_YSTRIPES + + + +**Return** + +- A geometry or an error. Use [tg_geom_error()](#group___geometry_parsing_1gae77b27ad34c2a215cc281647ab6dbc7e) after parsing to check for errors. + + +**See also** + +- [tg_parse_geobin_ix()](#group___geometry_parsing_1gac0450996bcd71cdde81af451dd0e9571) +- [tg_geom_error()](#group___geometry_parsing_1gae77b27ad34c2a215cc281647ab6dbc7e) +- [tg_geom_geobin()](#group___geometry_writing_1gab0f92319db4b8c61c62e0bd735d907b7) +- [https://github.com/tidwall/tg/blob/main/docs/GeoBIN.md](https://github.com/tidwall/tg/blob/main/docs/GeoBIN.md) +- [Geometry parsing](#group___geometry_parsing) + + + + +## tg_parse_geobin_ix() +```c +struct tg_geom *tg_parse_geobin_ix(const uint8_t *geobin, size_t len, enum tg_index ix); +``` +Parse GeoBIN binary using provided indexing option. + +**Parameters** + +- **geobin**: GeoBIN data +- **len**: Length of data +- **ix**: Indexing option, e.g. TG_NONE, TG_NATURAL, TG_YSTRIPES + + + +**Return** + +- A geometry or an error. Use [tg_geom_error()](#group___geometry_parsing_1gae77b27ad34c2a215cc281647ab6dbc7e) after parsing to check for errors. + + +**See also** + +- [tg_parse_geobin()](#group___geometry_parsing_1ga4792bf4f319356fff3fa539cbd44b196) + + + ## tg_parse() ```c @@ -2291,6 +2378,34 @@ if (tg_geom_error(geom)) { + +## tg_geobin_fullrect() +```c +int tg_geobin_fullrect(const uint8_t *geobin, size_t len, double min[4], double max[4]); +``` +Returns the minimum bounding rectangle of GeoBIN data. + +**Parameters** + +- **geobin**: GeoBIN data +- **len**: Length of data +- **min**: min values, must have room for 4 dimensions +- **max**: max values, must have room for 4 dimensions + + + +**Return** + +- number of dimensions, or zero if rect cannot be determined. + + +**See also** + +- [tg_geom_fullrect()](#group___geometry_accessors_1gac1a077f09e247c022f09e48392b80051) +- [tg_geom_rect()](#group___geometry_accessors_1ga68d67f900b847ae08e6515a620f4f657) + + + ## tg_geom_geojson() ```c @@ -2459,6 +2574,49 @@ if (len > sizeof(str)-1) { + +## tg_geom_geobin() +```c +size_t tg_geom_geobin(const struct tg_geom *geom, uint8_t *dst, size_t n); +``` +Writes a GeoBIN representation of a geometry. + +The content is stored in the buffer pointed by dst. + + + +**Parameters** + +- **geom**: Input geometry +- **dst**: Buffer where the resulting content is stored. +- **n**: Maximum number of bytes to be used in the buffer. + + + +**Return** + +- The number of characters needed to store the content into the buffer. If the returned length is greater than n, then only a parital copy occurred, for example: +```c +uint8_t buf[64]; +size_t len = tg_geom_geobin(geom, buf, sizeof(buf)); +if (len > sizeof(buf)) { + // ... write did not complete ... +} +``` + + + + +**See also** + +- [tg_geom_geojson()](#group___geometry_writing_1gae4669b8ee598a46e7f4c5a11e645fe8f) +- [tg_geom_wkt()](#group___geometry_writing_1ga047599bbe51886a6c3fbe35a6372d5d6) +- [tg_geom_wkb()](#group___geometry_writing_1gac0938697f9270d1924c6746af68dd693) +- [tg_geom_hex()](#group___geometry_writing_1ga269f61715302d48f144ffaf46e8c397b) +- [Geometry writing](#group___geometry_writing) + + + ## tg_geom_new_point_z() ```c diff --git a/docs/GEOBIN.md b/docs/GEOBIN.md new file mode 100644 index 0000000..4ce8301 --- /dev/null +++ b/docs/GEOBIN.md @@ -0,0 +1,137 @@ +# GeoBIN + +A simple binary format for GeoJSON. + +## Features + +- Geometry data uses standard WKB format +- Includes the geometry bounds (MBR) +- All JSON members are preserved +- Only little endian WKB allowed + +#### Goals + +To provide a simple lossless format for storing GeoJSON as binary, and to have +the MBR easily accessible for spatial indexing. + +It's a non-goal to have an emphasis on minimizing sizes, other that what is +already provided by converting JSON into WKB. In fact, GeoBIN will be a tiny +bit larger than WKB for non-point types. Roughly 34 bytes larger in most cases. + +## The Format + +The GeoBIN format consists of the Head byte, MBR Section, Extra JSON Section, +and the geometries WKB. In the case of a FeatureCollection, each Feature is +also represented as GeoBIN. + +### Head byte + +- 0x01: Point WKB (includes this byte) +- 0x02: Non-point Geometry +- 0x03: GeoJSON Feature +- 0x04: GeoJSON FeatureCollection + +### MBR Section + +All binary except 0x01 has an MBR that immediately follows the head byte. + +The MBR format starts with a single byte integer that is the number of +dimensions, followed by a series of coordinates that define the rectangle. + +#### 2D MBR + +- 1-byte integer: number of dimensions (2) +- 8-byte float[4]: `[xmin, ymin, xmax, ymax]` + +#### 3D MBR + +- 1-byte integer: number of dimensions (3) +- 8-byte float[6]: `[xmin, ymin, zmin, xmax, ymax, zmax]` + +#### 4D MBR + +- 1-byte integer: number of dimensions (4) +- 8-byte float[8]: `[xmin, ymin, zmin, mmin, xmax, ymax, zmax, mmax]` + +GeoBIN only supports 2, 3, or 4 dimensions; or zero to omit this section +altogether. + +### Extra JSON Section + +GeoJSON is allowed to have any number of extra json members, other than the +standard "type", "coordinates", "geometries", and "features". These extra json +members must be maintained in order to preserve the original GeoJSON object. + +Immediately after the MBR Section is the Extra JSON Section, which is a +null-terminated string. + +For example given this GeoJSON. + +```json +{ + "type": "Feature", + "id": 1934, + "geometry": { "type": "Point", "coordinates": [-112, 33] }, + "properties": { + "terrain": "desert" + } +} +``` + +The "id" and "properties" are the 'extra json', and should be extracted into +its own mini json document. + +```json +{"id":1934,"properties":{"terrain":"desert"}} +``` + +An empty string is just the '\0' null character indicating that there is no +extra json. + +It's required that when provided the extra json is valid json. + +### 0x01: Point WKB + +When the head byte is 0x01, the entire binary is a WKB Point. +The MBR must be cacluated from the WKB, which shouldn't be too difficult since +the WKB has its own head byte specifying dimensionality. + +### 0x02: Non-point Geometry + +When the head byte is 0x02, the binary will have an MBR, Extra JSON, and WKB. +Immediately after the Extra JSON Section is the raw geometry WKB. + +### 0x03: Feature + +Just like a 0x02, but indicates that the geometry is wrapped in a Feature +object. + +The following two are effictively the same thing, except that the first uses +0x02 and the second uses 0x03. + +```json +{ "type": "LineString", "coordinates": [[10, 10], [20, 20]] } +``` + +```json +{ + "type": "Feature", + "geometry": { "type": "LineString", "coordinates": [[10, 10], [20, 20]] } +} +``` + +### 0x04: FeatureCollection + +The 0x04 head byte means that the original GeoJSON is a FeatureCollection. +It's just like 0x02 except that after the Extra JSON Section, in lieu of raw +WKB, is a list of Features. + +This list starts with the count, followed by the features. Each feature is a +complete GeoBin object. + +- 4-byte integer: number of features +- N-features: series of features. Each a valid GeoBIN object. + +--- + +GeoBIN support is provided by [TG](https://github.com/tidwall/tg). diff --git a/tests/test_geobin.c b/tests/test_geobin.c new file mode 100644 index 0000000..eaee409 --- /dev/null +++ b/tests/test_geobin.c @@ -0,0 +1,173 @@ +#include "tests.h" + +void perfect_match(const char *geojson_or_wkt, const char *expect, int dims, double min[4], double max[4]) { + struct tg_geom *g1; + if (geojson_or_wkt[0] == '{') { + g1 = tg_parse_geojson(geojson_or_wkt); + } else { + g1 = tg_parse_wkt(geojson_or_wkt); + } + if (tg_geom_error(g1)) { + printf("%s\n", tg_geom_error(g1)); + assert(0); + } + + double gmin[4], gmax[4]; + int gdims = tg_geom_fullrect(g1, gmin, gmax); + cmpfullrect(gdims, gmin, gmax, dims, min, max); + uint8_t geobin[100000]; + size_t geobinlen = tg_geom_geobin(g1, geobin, sizeof(geobin)); + assert(geobinlen <= sizeof(geobin)); + struct tg_geom *g2 = tg_parse_geobin(geobin, geobinlen); + if (tg_geom_error(g2)) { + printf("%s\n", tg_geom_error(g2)); + assert(0); + } + + + // printf("0x%02X\n", geobin[0]); + gdims = tg_geobin_fullrect(geobin, geobinlen, gmin, gmax); + + cmpfullrect(gdims, gmin, gmax, dims, min, max); + + char geojson_or_wkt2[100000]; + size_t geojson_or_wkt2len; + if (geojson_or_wkt[0] == '{') { + geojson_or_wkt2len = tg_geom_geojson(g2, geojson_or_wkt2, sizeof(geojson_or_wkt2)); + } else { + geojson_or_wkt2len = tg_geom_wkt(g2, geojson_or_wkt2, sizeof(geojson_or_wkt2)); + } + assert(geojson_or_wkt2len < sizeof(geojson_or_wkt2)); + + if (!expect) { + expect = geojson_or_wkt; + } + if (strcmp(geojson_or_wkt2, expect) != 0) { + printf("expected:\n\t%s\ngot\n\t%s\n", expect, geojson_or_wkt2); + assert(0); + } + + tg_geom_free(g1); + tg_geom_free(g2); +} + +void test_geobin_basic_syntax(void) { + assert(tg_geom_geobin(0, 0, 0) == 0); + perfect_match("{\"type\":\"Point\",\"coordinates\":[55,44]}", 0, 2, (double[]){55, 44}, (double[]){55, 44}); + perfect_match("{\"type\":\"Point\",\"coordinates\":[55,44,33]}", 0, 3, (double[]){55, 44, 33}, (double[]){55, 44, 33}); + perfect_match("{\"type\":\"Point\",\"coordinates\":[55,44,33,22]}", 0, 4, (double[]){55, 44, 33, 22}, (double[]){55, 44, 33, 22}); + perfect_match("{\"type\":\"Point\",\"coordinates\":[55,44],\"a\":{\"b\":\"c\"}}", 0, 2, (double[]){55, 44}, (double[]){55, 44}); + perfect_match("{\"type\":\"Point\",\"coordinates\":[55,44,33],\"a\":{\"b\":\"c\"}}", 0, 3, (double[]){55, 44, 33}, (double[]){55, 44, 33}); + perfect_match("{\"type\":\"Point\",\"coordinates\":[55,44,33,22],\"a\":{\"b\":\"c\"}}", 0, 4, (double[]){55, 44, 33, 22}, (double[]){55, 44, 33, 22}); + perfect_match("POINT(55 44)", 0, 2, (double[]){55, 44}, (double[]){55, 44}); + perfect_match("POINT(55 44 33)", 0, 3, (double[]){55, 44, 33}, (double[]){55, 44, 33}); + perfect_match("POINT(55 44 33 22)", 0, 4, (double[]){55, 44, 33, 22}, (double[]){55, 44, 33, 22}); + perfect_match("POINT(55 44)", 0, 2, (double[]){55, 44}, (double[]){55, 44}); + perfect_match("POINTZ(55 44 33)", "POINT(55 44 33)", 3, (double[]){55, 44, 33}, (double[]){55, 44, 33}); + perfect_match("POINTM(55 44 33)", "POINT M(55 44 33)", 3, (double[]){55, 44, 33}, (double[]){55, 44, 33}); + perfect_match("POINT(55 44 33 22)", 0, 4, (double[]){55, 44, 33, 22}, (double[]){55, 44, 33, 22}); + perfect_match("{\"type\":\"LineString\",\"coordinates\":[[55,44],[33,22]]}", 0, 2, (double[]){33, 22}, (double[]){55, 44}); + perfect_match("{\"type\":\"LineString\",\"coordinates\":[[55,44,11],[33,22,66]]}", 0, 3, (double[]){33, 22,11}, (double[]){55, 44,66}); + perfect_match("{\"type\":\"LineString\",\"coordinates\":[[55,44,11,99],[33,22,66,77]]}", 0, 4, (double[]){33, 22, 11, 77}, (double[]){55, 44, 66, 99}); + perfect_match("{\"type\":\"Polygon\",\"coordinates\":[[[33,44],[55,66],[33,44]]]}", 0, 2, (double[]){33, 44}, (double[]){55, 66}); + perfect_match("{\"type\":\"Polygon\",\"coordinates\":[[[33,44],[55,66],[44,66],[33,44]],[[42,61],[46,61],[46,64],[42,61]]]}", 0, 2, (double[]){33, 44}, (double[]){55, 66}); + perfect_match( + "{\"type\":\"Feature\",\"geometry\":{\"type\":\"Polygon\",\"coordinates\":[[[33,44],[55,66],[44,66],[33,44]],[[42,61],[46,61],[46,64],[42,61]]]}}", + "{\"type\":\"Feature\",\"geometry\":{\"type\":\"Polygon\",\"coordinates\":[[[33,44],[55,66],[44,66],[33,44]],[[42,61],[46,61],[46,64],[42,61]]]},\"properties\":{}}", + 2, (double[]){33, 44}, (double[]){55, 66}); + perfect_match( + "{\"type\":\"Feature\",\"id\":123,\"geometry\":{\"type\":\"Polygon\",\"coordinates\":[[[33,44],[55,66],[44,66],[33,44]],[[42,61],[46,61],[46,64],[42,61]]]},\"properties\":{\"a\":\"b\"}}", + 0, 2, (double[]){33, 44}, (double[]){55, 66}); + perfect_match( + "{\"type\":\"FeatureCollection\",\"features\":[" + "{\"type\":\"Feature\",\"id\":123,\"geometry\":{\"type\":\"Polygon\",\"coordinates\":[[[33,44],[55,66],[44,66],[33,44]],[[42,61],[46,61],[46,64],[42,61]]]},\"properties\":{\"a\":\"b\"}}," + "{\"type\":\"Feature\",\"id\":123,\"geometry\":{\"type\":\"Polygon\",\"coordinates\":[[[33,44],[55,66],[44,66],[33,44]],[[42,61],[46,61],[46,64],[42,61]]]},\"properties\":{\"a\":\"b\"}}" + "]}", + 0, 2, (double[]){33, 44}, (double[]){55, 66}); +} + +void test_geobin_fail() { + assert(tg_geom_error(gc_geom(tg_parse_geobin(0, 0)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x00},1)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x05},1)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x02}, 1)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x02,2}, 2)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x02,1}, 2)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x02,5}, 2)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x02,0}, 2)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x02,0,0,0x1}, 4)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x04,0,0,0x1,0x0,0x0}, 6)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x04,0,0,0x1,0x0,0x0,0x0}, 7)))); + assert(tg_geom_error(gc_geom(tg_parse((uint8_t[]){0x04,0,0,0x1,0x0,0x0,0x0,0x2,0x0,0x0}, 10)))); + // invalid extra json. just the '{' character + assert(tg_geom_error(gc_geom(tg_parse_hex("02007B000101000000000000000000F87F000000000000F87F")))); +} + +void test_geobin_chaos() { + uint8_t geobin[5000]; + size_t len; + struct tg_geom *g; + while (1) { + g = tg_parse_geojson( + "{\"type\":\"FeatureCollection\",\"features\":[" + "{\"type\":\"Feature\",\"id\":123,\"geometry\":{\"type\":\"Polygon\",\"coordinates\":[[[33,44],[55,66],[44,66],[33,44]],[[42,61],[46,61],[46,64],[42,61]]]},\"properties\":{\"a\":\"b\"}}," + "{\"type\":\"Feature\",\"id\":123,\"geometry\":{\"type\":\"Polygon\",\"coordinates\":[[[33,44],[55,66],[44,66],[33,44]],[[42,61],[46,61],[46,64],[42,61]]]},\"properties\":{\"a\":\"b\"}}" + "]}"); + len = tg_geom_geobin(g, geobin, sizeof(geobin)); + if (g) break; + } + double secs = 3.0; + double start = now(); + while (now()-start < secs) { + struct tg_geom *g2 = tg_parse_geobin(geobin, len); + tg_geom_free(g2); + } + tg_geom_free(g); +} + +char *make_deep_geobin_hex(int depth) { + char gcname[] = "04000001000000"; // FeatureCollection() one child + char ptname[] = "0101000000000000000000F87F000000000000F87F"; // POINT EMPTY + char *hex = malloc((strlen(gcname)*depth)+strlen(ptname)+1); + assert(hex); + int len = 0; + // memcpy(hex+len, gbin, strlen(gbin)); + // len += strlen(gbin); + for (int i = 0;i < depth; i++) { + memcpy(hex+len, gcname, strlen(gcname)); + len += strlen(gcname); + } + memcpy(hex+len, ptname, strlen(ptname)); + len += strlen(ptname); + hex[len] = '\0'; + return hex; +} + +void test_geobin_max_depth() { + // nested geometry collections + char *hex; + struct tg_geom *geom; + + int depths[] = { 1, 100, 1000, 1023, 1024, 1025, 2000 }; + + for (size_t i = 0; i < sizeof(depths)/sizeof(int); i++) { + hex = make_deep_geobin_hex(depths[i]); + geom = tg_parse_hex(hex); + if (depths[i] > 1024) { + assert(tg_geom_error(geom)); + } else { + assert(!tg_geom_error(geom)); + } + tg_geom_free(geom); + free(hex); + } +} + + +int main(int argc, char **argv) { + do_test(test_geobin_basic_syntax); + do_test(test_geobin_fail); + do_test(test_geobin_max_depth); + do_chaos_test(test_geobin_chaos); + return 0; +} diff --git a/tests/test_rect.c b/tests/test_rect.c index 2e2e016..33d57d9 100644 --- a/tests/test_rect.c +++ b/tests/test_rect.c @@ -159,6 +159,68 @@ void test_rect_distance() { assert(tg_rect_distance_rect(R(11,10,11,10), R(0,0,10,10)) == 1); } + +void fullrectmatch(const char *wkt, int dims, double min[4], double max[4]) { + struct tg_geom *g = tg_parse_wkt(wkt); + if (tg_geom_error(g)) { + printf("%s\n", tg_geom_error(g)); + assert(0); + } + + double min2[4], max2[4]; + int dims2 = tg_geom_fullrect(g, min2, max2); + tg_geom_free(g); + cmpfullrect(dims, min, max, dims2, min2, max2); +} + +void test_rect_fullrect() { + assert(tg_geom_fullrect(0, 0, 0) == 0); + fullrectmatch("POINT(10 20)", 2, (double[]){10, 20}, (double[]){10, 20}); + fullrectmatch("POINT(10 20 30)", 3, (double[]){10, 20, 30}, (double[]){10, 20, 30}); + fullrectmatch("POINT(10 20 30 40)", 4, (double[]){10, 20, 30, 40}, (double[]){10, 20, 30, 40}); + fullrectmatch("LINESTRING(10 20,30 40)", 2,(double[]){10,20},(double[]){30,40}); + fullrectmatch("LINESTRING(10 20 30,40 50 60)", 3,(double[]){10,20,30},(double[]){40,50,60}); + fullrectmatch("LINESTRING(40 50 60,10 20 30)", 3,(double[]){10,20,30},(double[]){40,50,60}); + fullrectmatch("LINESTRING(10 20 30 40,50 60 70 80)", 4,(double[]){10,20,30,40},(double[]){50,60,70,80}); + fullrectmatch("LINESTRING(50 60 70 80,10 20 30 40)", 4,(double[]){10,20,30,40},(double[]){50,60,70,80}); + fullrectmatch("GEOMETRYCOLLECTION(" + "POINT(10 20)" + ")", 2, (double[]){10, 20}, (double[]){10, 20}); + fullrectmatch("GEOMETRYCOLLECTION(" + "POINT(10 20)," + "POINT(30 40)" + ")", 2, (double[]){10, 20}, (double[]){30, 40}); + fullrectmatch("GEOMETRYCOLLECTION(" + "POINT(10 20 30)," + "POINT(40 50 60)" + ")", 3, (double[]){10, 20, 30}, (double[]){40, 50, 60}); + fullrectmatch("GEOMETRYCOLLECTION(" + "POINT(10 20 30 40)," + "POINT(50 60 70 80)" + ")", 4, (double[]){10, 20, 30, 40}, (double[]){50, 60, 70, 80}); + fullrectmatch("GEOMETRYCOLLECTION(" + "POINT(40 50 60)," + "POINT(10 20 30)" + ")", 3, (double[]){10, 20, 30}, (double[]){40, 50, 60}); + fullrectmatch("GEOMETRYCOLLECTION(" + "POINT(50 60 70 80)," + "POINT(10 20 30 40)" + ")", 4, (double[]){10, 20, 30, 40}, (double[]){50, 60, 70, 80}); + fullrectmatch("GEOMETRYCOLLECTION(" + "POINT(10 20)," + "POINT(30 40 50)" + ")", 3, (double[]){10, 20, 50}, (double[]){30, 40, 50}); + fullrectmatch("GEOMETRYCOLLECTION(" + "POINT(10 20)," + "POINT(30 40 50 60)" + ")", 4, (double[]){10, 20, 50, 60}, (double[]){30, 40, 50, 60}); + fullrectmatch("GEOMETRYCOLLECTION(" + "POINT(10 20)," + "POINT(30 40 50 60)," + "POINT(70 80 90)" + ")", 4, (double[]){10, 20, 50, 60}, (double[]){70, 80, 90, 60}); +} + int main(int argc, char **argv) { do_test(test_rect_center); do_test(test_rect_move); @@ -176,4 +238,5 @@ int main(int argc, char **argv) { do_test(test_rect_covers_poly); do_test(test_rect_intersects_poly); do_test(test_rect_distance); + do_test(test_rect_fullrect); } diff --git a/tests/tests.h b/tests/tests.h index d77aab0..3bf7f8b 100644 --- a/tests/tests.h +++ b/tests/tests.h @@ -658,6 +658,34 @@ char *read_file(const char *path, size_t *size) { return data; } +void cmpfullrect(int dims, double min[4], double max[4], int dims2, double min2[4], double max2[4]) { + if (dims != dims2) { + printf("expected\n\tdims( %d )\ngot\n\tdims( %d )\n", dims, dims2); + assert(0); + } + for (int i = 0; i < dims; i++) { + if (min[i] != min2[i] || max[i] != max2[i]) { + printf("expected\n\tmin( "); + for (int i = 0; i < dims; i++) { + printf("%.0f ", min2[i]); + } + printf(") max( "); + for (int i = 0; i < dims; i++) { + printf("%.0f ", max2[i]); + } + printf(")\ngot\n\tmin( "); + for (int i = 0; i < dims; i++) { + printf("%.0f ", min[i]); + } + printf(") max( "); + for (int i = 0; i < dims; i++) { + printf("%.0f ", max[i]); + } + printf(")\n"); + assert(0); + } + } +} static double *make_polygon(int n, double x, double y, double r, double c, double s) { s = s < 0.0 ? 0.0 : s > 1.0 ? 1.0 : s; diff --git a/tg.c b/tg.c index f4d017f..11bf6e6 100644 --- a/tg.c +++ b/tg.c @@ -8298,7 +8298,7 @@ static void write_doublele(struct writer *wr, double x) { } static void write_string(struct writer *wr, const char *str) { - char *p = (char*)str; + const char *p = str; while (*p) write_char(wr, *(p++)); } @@ -11105,15 +11105,6 @@ static int parse_wkt_posns(enum base base, int dims, int depth, const char *wkt, if (xparens) { if (i == len || wkt[i] != '(') { // err: expected '(' - // if (wkt[i] == 'e' || wkt[i] == 'E') { - // (wkt[i+1] == 'm' || wkt[i+1] == 'M') - // (wkt[i+2] == 'p' || wkt[i+2] == 'P') - // (wkt[i+3] == 't' || wkt[i+3] == 'T') - // (wkt[i+4] == 'y' || wkt[i+4] == 'Y') - // if (i+5 - // // EMPTY ? - // } - // printf(">>> %c\n", wkt[i]); *err = wkt_invalid_err("expected '('"); return -1; } @@ -13439,6 +13430,9 @@ size_t tg_geom_hex(const struct tg_geom *geom, char *dst, size_t n) { return count*2; } +static size_t parse_geobin(const uint8_t *geobin, size_t len, size_t i, + size_t depth, enum tg_index ix, struct tg_geom **g); + static struct tg_geom *parse_hex(const char *hex, size_t len, enum tg_index ix) { const uint8_t _ = 0; @@ -13468,7 +13462,14 @@ static struct tg_geom *parse_hex(const char *hex, size_t len, enum tg_index ix) j++; } struct tg_geom *geom; - parse_wkb(dst, len/2, 0, 0, ix, &geom); + len /= 2; + size_t n; + if (len > 0 && dst[0] >= 0x2 && dst[0] <= 0x4) { + n = parse_geobin(dst, len, 0, 0, ix, &geom); + } else { + n = parse_wkb(dst, len, 0, 0, ix, &geom); + } + (void)n; if (must_free) tg_free(dst); return geom; invalid: @@ -13476,7 +13477,8 @@ static struct tg_geom *parse_hex(const char *hex, size_t len, enum tg_index ix) return make_parse_error(wkb_invalid_err()); } -/// Parse hex encoded Well-known binary (WKB) using provided indexing option. +/// Parse hex encoded Well-known binary (WKB) or GeoBIN using provided indexing +/// option. /// @param hex Hex data /// @param len Length of data /// @param ix Indexing option, e.g. TG_NONE, TG_NATURAL, TG_YSTRIPES @@ -13498,7 +13500,8 @@ struct tg_geom *tg_parse_hexn_ix(const char *hex, size_t len, return geom; } -/// Parse hex encoded Well-known binary (WKB) using provided indexing option. +/// Parse hex encoded Well-known binary (WKB) or GeoBIN using provided indexing +/// option. /// @param hex Hex string. Must be null-terminated /// @param ix Indexing option, e.g. TG_NONE, TG_NATURAL, TG_YSTRIPES /// @returns A geometry or an error. Use tg_geom_error() after parsing to check @@ -13510,7 +13513,8 @@ struct tg_geom *tg_parse_hex_ix(const char *hex, enum tg_index ix) { return tg_parse_hexn_ix(hex, hex?strlen(hex):0, ix); } -/// Parse hex encoded Well-known binary (WKB) with an included data length. +/// Parse hex encoded Well-known binary (WKB) or GeoBIN with an included data +/// length. /// @param hex Hex data /// @param len Length of data /// @returns A geometry or an error. Use tg_geom_error() after parsing to check @@ -13521,7 +13525,7 @@ struct tg_geom *tg_parse_hexn(const char *hex, size_t len) { return tg_parse_hexn_ix(hex, len, TG_DEFAULT); } -/// Parse hex encoded Well-known binary (WKB). +/// Parse hex encoded Well-known binary (WKB) or GeoBIN. /// @param hex A hex string. Must be null-terminated /// @returns A geometry or an error. Use tg_geom_error() after parsing to check /// for errors. @@ -14352,7 +14356,10 @@ struct tg_geom *tg_parse_ix(const void *data, size_t len, enum tg_index ix) { } goto wkt; } - goto wkb; + if (src[0] == 0x00 || src[0] == 0x01) { + goto wkb; + } + goto geobin; geojson: return tg_parse_geojsonn_ix(src, len, ix); wkt: @@ -14361,6 +14368,8 @@ struct tg_geom *tg_parse_ix(const void *data, size_t len, enum tg_index ix) { return tg_parse_hexn_ix(src, len, ix); wkb: return tg_parse_wkb_ix((uint8_t*)src, len, ix); +geobin: + return tg_parse_geobin_ix((uint8_t*)src, len, ix); } /// Utility for returning an error message wrapped in a geometry. @@ -14377,3 +14386,426 @@ void tg_geom_setnoheap(struct tg_geom *geom) { geom->head.rc = 0; geom->head.noheap = 1; } + +/// Parse GeoBIN binary using provided indexing option. +/// @param geobin GeoBIN data +/// @param len Length of data +/// @param ix Indexing option, e.g. TG_NONE, TG_NATURAL, TG_YSTRIPES +/// @returns A geometry or an error. Use tg_geom_error() after parsing to check +/// for errors. +/// @see tg_parse_geobin_ix() +/// @see tg_geom_error() +/// @see tg_geom_geobin() +/// @see https://github.com/tidwall/tg/blob/main/docs/GeoBIN.md +/// @see GeometryParsing +struct tg_geom *tg_parse_geobin(const uint8_t *geobin, size_t len) { + return tg_parse_geobin_ix(geobin, len, 0); +} + +static size_t parse_geobin(const uint8_t *geobin, size_t len, size_t i, + size_t depth, enum tg_index ix, struct tg_geom **g) +{ + if (i == len) goto invalid; + if (depth > MAXDEPTH) goto invalid; + int head = geobin[i]; + if (head == 0x01) { + return parse_wkb(geobin, len, i, depth, ix, g); + } + i++; + if (head < 0x02 || head > 0x04) { + goto invalid; + } + if (i == len) { + goto invalid; + } + int dims = geobin[i++]; + if (dims && (dims < 2 || dims > 4)) { + goto invalid; + } + if (dims) { + i += 8*dims*2; + if (i >= len) { + goto invalid; + } + } + size_t xjsonlen = 0; + const char *xjson = (const char*)(geobin+i); + for (; i < len; i++) { + if (geobin[i] == '\0') { + i++; + break; + } + xjsonlen++; + } + if (i == len) { + goto invalid; + } + if (xjsonlen > 0 && !json_validn(xjson, xjsonlen)) { + goto invalid; + } + struct tg_geom *geom; + if (head == 0x04) { + // FeatureCollection + if (i+4 > len) { + goto invalid; + } + uint32_t nfeats; + memcpy(&nfeats, geobin+i, 4); + i += 4; + struct tg_geom **feats = tg_malloc(nfeats*sizeof(struct tg_geom*)); + if (!feats) { + return 0; + } + struct tg_geom *feat = 0; + uint32_t j = 0; + for (; j < nfeats; j++) { + i = parse_geobin(geobin, len, i, depth+1, ix, &feat); + if (i == PARSE_FAIL) { + break; + } + feats[j] = feat; + } + if (j == nfeats) { + geom = tg_geom_new_geometrycollection((void*)feats, nfeats); + } + for (uint32_t k = 0; k < j; k++) { + tg_geom_free(feats[k]); + } + tg_free(feats); + if (j < nfeats) { + *g = feat; // return the last failed feature + return PARSE_FAIL; + } + if (!geom) { + *g = 0; + return PARSE_FAIL; + } + geom->head.flags |= IS_FEATURE_COL; + } else { + i = parse_wkb(geobin, len, i, depth, ix, &geom); + } + if (i == PARSE_FAIL || !geom) { + *g = geom; + return PARSE_FAIL; + } + if ((xjsonlen > 0 || head == 0x03) && geom->head.base != BASE_GEOM) { + // Wrap base in tg_geom + struct tg_geom *g2 = geom_new(geom->head.type); + if (!g2) { + tg_geom_free(geom); + *g = 0; + return PARSE_FAIL; + } + if (geom->head.base == BASE_POINT) { + g2->point = ((struct boxed_point*)geom)->point; + boxed_point_free((struct boxed_point*)geom); + } else { + g2->line = (struct tg_line*)geom; + } + geom = g2; + } + if (head == 0x03) { + geom->head.flags |= IS_FEATURE; + } + if (xjsonlen > 0) { + geom->xjson = tg_malloc(xjsonlen+1); + if (!geom->xjson) { + tg_geom_free(geom); + *g = 0; + return PARSE_FAIL; + } + memcpy(geom->xjson, xjson, xjsonlen+1); + } + *g = geom; + return i; +invalid: + *g = make_parse_error("invalid binary"); + return PARSE_FAIL; + +} + +/// Parse GeoBIN binary using provided indexing option. +/// @param geobin GeoBIN data +/// @param len Length of data +/// @param ix Indexing option, e.g. TG_NONE, TG_NATURAL, TG_YSTRIPES +/// @returns A geometry or an error. Use tg_geom_error() after parsing to check +/// for errors. +/// @see tg_parse_geobin() +struct tg_geom *tg_parse_geobin_ix(const uint8_t *geobin, size_t len, + enum tg_index ix) +{ + struct tg_geom *geom = NULL; + parse_geobin(geobin, len, 0, 0, ix, &geom); + if (!geom) return NULL; + if ((geom->head.flags&IS_ERROR) == IS_ERROR) { + struct tg_geom *gerr = make_parse_error("ParseError: %s", geom->error); + tg_geom_free(geom); + return gerr; + } + return geom; +} + +static void write_geom_geobin(const struct tg_geom *geom, struct writer *wr); + +static void write_base_geom_geobin(const struct tg_geom *geom, + struct writer *wr) +{ + // extra json section + const char *xjson = tg_geom_extra_json(geom); + + // write head byte + if ((geom->head.flags&IS_FEATURE_COL) == IS_FEATURE_COL) { + write_byte(wr, 0x04); + } else if ((geom->head.flags&IS_FEATURE) == IS_FEATURE) { + write_byte(wr, 0x03); + } else if (geom->head.type == TG_POINT && !xjson) { + write_geom_point_wkb(geom, wr); + return; + } else { + write_byte(wr, 0x02); + } + + // mbr section + double min[4], max[4]; + int dims = tg_geom_fullrect(geom, min, max); + write_byte(wr, dims); + for (int i = 0; i < dims; i++) { + write_doublele(wr, min[i]); + } + for (int i = 0; i < dims; i++) { + write_doublele(wr, max[i]); + } + + if (xjson) { + write_string(wr, xjson); + } + write_byte(wr, 0); + + if ((geom->head.flags&IS_FEATURE_COL) == IS_FEATURE_COL) { + // write feature collection + int ngeoms = tg_geom_num_geometries(geom); + write_uint32le(wr, (uint32_t)ngeoms); + for (int i = 0; i < ngeoms; i++) { + const struct tg_geom *g2 = tg_geom_geometry_at(geom, i); + write_geom_geobin(g2, wr); + } + } else { + // write wkb + write_geom_wkb(geom, wr); + } +} + +static void write_point_geobin(struct boxed_point *point, struct writer *wr) { + write_point_wkb(point, wr); +} + +static void write_geobin_rect(struct writer *wr, struct tg_rect rect) { + write_byte(wr, 2); // dims + write_doublele(wr, rect.min.x); + write_doublele(wr, rect.min.y); + write_doublele(wr, rect.max.x); + write_doublele(wr, rect.max.y); +} + +static void write_line_geobin(struct tg_line *line, struct writer *wr) { + write_byte(wr, 0x02); + write_geobin_rect(wr, ((struct tg_ring*)line)->rect); + write_byte(wr, 0); + write_line_wkb(line, wr); +} + +static void write_ring_geobin(struct tg_ring *ring, struct writer *wr) { + write_byte(wr, 0x02); + write_geobin_rect(wr, ring->rect); + write_byte(wr, 0); + write_ring_wkb(ring, wr); +} + +static void write_poly_geobin(struct tg_poly *poly, struct writer *wr) { + write_byte(wr, 0x02); + write_geobin_rect(wr, poly->exterior->rect); + write_byte(wr, 0); + write_poly_wkb(poly, wr); +} + +static void write_geom_geobin(const struct tg_geom *geom, struct writer *wr) { + if ((geom->head.flags&IS_FEATURE) == IS_FEATURE) { + goto base_geom; + } + switch (geom->head.base) { + case BASE_GEOM: + base_geom: + write_base_geom_geobin(geom, wr); + break; + case BASE_POINT: + write_point_geobin((struct boxed_point*)geom, wr); + break; + case BASE_LINE: + write_line_geobin((struct tg_line*)geom, wr); + break; + case BASE_RING: + write_ring_geobin((struct tg_ring*)geom, wr); + break; + case BASE_POLY: + write_poly_geobin((struct tg_poly*)geom, wr); + break; + } +} + +/// Writes a GeoBIN representation of a geometry. +/// +/// The content is stored in the buffer pointed by dst. +/// +/// @param geom Input geometry +/// @param dst Buffer where the resulting content is stored. +/// @param n Maximum number of bytes to be used in the buffer. +/// @return The number of characters needed to store the content into the +/// buffer. +/// If the returned length is greater than n, then only a parital copy +/// occurred, for example: +/// +/// ``` +/// uint8_t buf[64]; +/// size_t len = tg_geom_geobin(geom, buf, sizeof(buf)); +/// if (len > sizeof(buf)) { +/// // ... write did not complete ... +/// } +/// ``` +/// +/// @see tg_geom_geojson() +/// @see tg_geom_wkt() +/// @see tg_geom_wkb() +/// @see tg_geom_hex() +/// @see GeometryWriting +size_t tg_geom_geobin(const struct tg_geom *geom, uint8_t *dst, size_t n) { + if (!geom) return 0; + struct writer wr = { .dst = dst, .n = n }; + write_geom_geobin(geom, &wr); + return wr.count; +} + +/// Returns the minimum bounding rectangle of a geometry on all dimensions. +/// @param geom Input geometry +/// @param min min values, must have room for 4 dimensions +/// @param max max values, must have room for 4 dimensions +/// @return number of dimensions, or zero if invalid geom. +/// @see tg_geom_rect() +int tg_geom_fullrect(const struct tg_geom *geom, double min[4], double max[4]) { + if (!geom) { + return 0; + } + struct tg_rect rect = tg_geom_rect(geom); + min[0] = rect.min.x; + min[1] = rect.min.y; + min[2] = 0; + min[3] = 0; + max[0] = rect.max.x; + max[1] = rect.max.y; + max[2] = 0; + max[3] = 0; + int dims = 2; + if (geom->head.base == BASE_GEOM) { + if (geom->head.type == TG_POINT) { + // Point + if ((geom->head.flags&HAS_Z) == HAS_Z) { + min[dims] = geom->z; + max[dims] = geom->z; + dims++; + } + if ((geom->head.flags&HAS_M) == HAS_M) { + min[dims] = geom->m; + max[dims] = geom->m; + dims++; + } + } else if (geom->head.type == TG_GEOMETRYCOLLECTION && geom->multi) { + // GeometryCollection. Expand all child geometries + struct tg_geom **geoms = geom->multi->geoms; + int ngeoms = geom->multi->ngeoms; + double gmin[4], gmax[4]; + for (int i = 0; i < ngeoms; i++) { + int gdims = tg_geom_fullrect(geoms[i], gmin, gmax); + if (gdims >= 3) { + if (dims == 2) { + min[2] = gmin[2]; + max[2] = gmax[2]; + dims++; + } else { + min[2] = fmin0(min[2], gmin[2]); + max[2] = fmax0(max[2], gmax[2]); + } + } + if (gdims >= 4) { + if (dims == 3) { + min[3] = gmin[3]; + max[3] = gmax[3]; + dims++; + } else { + min[3] = fmin0(min[3], gmin[3]); + max[3] = fmax0(max[3], gmax[3]); + } + } + } + } else { + // Other geometries + if ((geom->head.flags&HAS_Z) == HAS_Z) dims++; + if ((geom->head.flags&HAS_M) == HAS_M) dims++; + if (dims == 3 && geom->ncoords > 0) { + min[2] = geom->coords[0]; + max[2] = geom->coords[0]; + for (int i = 1; i < geom->ncoords; i++) { + min[2] = fmin0(min[2], geom->coords[i]); + max[2] = fmax0(max[2], geom->coords[i]); + } + } else if (dims == 4 && geom->ncoords > 1) { + min[2] = geom->coords[0]; + min[3] = geom->coords[1]; + max[2] = geom->coords[0]; + max[3] = geom->coords[1]; + for (int i = 2; i < geom->ncoords-1; i+=2) { + min[2] = fmin0(min[2], geom->coords[i]); + min[3] = fmin0(min[3], geom->coords[i+1]); + max[2] = fmax0(max[2], geom->coords[i]); + max[3] = fmax0(max[3], geom->coords[i+1]); + } + } + } + } + return dims; +} + +/// Returns the minimum bounding rectangle of GeoBIN data. +/// @param geobin GeoBIN data +/// @param len Length of data +/// @param min min values, must have room for 4 dimensions +/// @param max max values, must have room for 4 dimensions +/// @return number of dimensions, or zero if rect cannot be determined. +/// @see tg_geom_fullrect() +/// @see tg_geom_rect() +int tg_geobin_fullrect(const uint8_t *geobin, size_t len, double min[4], + double max[4]) +{ + size_t dims = 0; + if (geobin && len > 2 && geobin[0] >= 0x01 && geobin[0] <= 0x04) { + if (geobin[0] == 0x01 && len >= 5) { + // Read Point + uint32_t type; + memcpy(&type, geobin+1, 4); + switch (type) { + case 1: dims = 2; break; + case 1001: dims = 3; break; + case 2001: dims = 3; break; + case 3001: dims = 4; break; + } + if (dims > 0 && len >= 5+8*dims) { + memcpy(min, geobin+5, 8*dims); + memcpy(max, geobin+5, 8*dims); + } + } else if (geobin[0] != 0x01 && len >= 2+8*geobin[1]*2){ + // Read MBR + dims = geobin[1]; + memcpy(min, geobin+2, 8*dims); + memcpy(max, geobin+2+8*dims, 8*dims); + } + } + return dims; +} diff --git a/tg.h b/tg.h index 175c991..d865384 100644 --- a/tg.h +++ b/tg.h @@ -122,6 +122,7 @@ size_t tg_geom_memsize(const struct tg_geom *geom); void tg_geom_search(const struct tg_geom *geom, struct tg_rect rect, bool (*iter)(const struct tg_geom *geom, int index, void *udata), void *udata); +int tg_geom_fullrect(const struct tg_geom *geom, double min[4], double max[4]); /// @} /// @defgroup GeometryPredicates Geometry predicates @@ -157,9 +158,13 @@ struct tg_geom *tg_parse_hex(const char *hex); struct tg_geom *tg_parse_hexn(const char *hex, size_t len); struct tg_geom *tg_parse_hex_ix(const char *hex, enum tg_index ix); struct tg_geom *tg_parse_hexn_ix(const char *hex, size_t len, enum tg_index ix); +struct tg_geom *tg_parse_geobin(const uint8_t *geobin, size_t len); +struct tg_geom *tg_parse_geobin_ix(const uint8_t *geobin, size_t len,enum tg_index ix); struct tg_geom *tg_parse(const void *data, size_t len); struct tg_geom *tg_parse_ix(const void *data, size_t len, enum tg_index ix); const char *tg_geom_error(const struct tg_geom *geom); +int tg_geobin_fullrect(const uint8_t *geobin, size_t len, double min[4], double max[4]); + /// @} /// @defgroup GeometryWriting Geometry writing @@ -169,6 +174,7 @@ size_t tg_geom_geojson(const struct tg_geom *geom, char *dst, size_t n); size_t tg_geom_wkt(const struct tg_geom *geom, char *dst, size_t n); size_t tg_geom_wkb(const struct tg_geom *geom, uint8_t *dst, size_t n); size_t tg_geom_hex(const struct tg_geom *geom, char *dst, size_t n); +size_t tg_geom_geobin(const struct tg_geom *geom, uint8_t *dst, size_t n); /// @} /// @defgroup GeometryConstructorsEx Geometry with alternative dimensions