geoid – Efficient storage and queriying of geographic data
(require geoid) | package: geoid |
The geoid library allows representing latidude/longitude coordinates using 64-bit integers. The integer values are ordered to preseve locality, which means that these IDs can be stored in databases and used to query and efficiently retrieve locations which are in a geographic region.
This library is inspired by the S2 Geometry library, in the sense that the integer IDs are determined in the same way. Until this libary has its own introduction section, you can read S2 Cell Hierarchy to understand how geoids work. What this library calls geoids, are called cells by the S2 libary.
Also, while inspired by the S2 library, this Racket module is an independent implementation that does not generate compatible IDs and does not aim to have the same API and functionality as that library.
A geoid is a 64 bit integer between first-valid-geoid and last-valid-geoid. In particular, 0 is not a valid geoid (they actualy start at 1)
Geoids which are close together represent geographic locations which are close together, but the reverse is not true: there are geographic locations which are close together, but their geoids are very different
The projection method, while not uniform accross the globe, it only has a small distortion and there are no singularities anywhere, including at the poles
Geoids split the earth surface at different levels: the highest level is level 30, where the earch surface is split into 6 faces. At each lower level, the geoids are split into four, for example at level 29, each of the 6 faces are split into four, producing 24 total geois. The subdivision goes on until level 0 is reached. At level 0, each geoid represents a patch of earth of approximately 0.7 cm2.
1 Converting to and from Latitude/Longitude
procedure
(lat-lng->geoid lat lng [#:level level]) → exact-integer?
lat : real? lng : real? level : exact-integer? = 0
(geoid->lat-lng geoid) →
real? real? geoid : exact-integer?
geoid->lat-lng will return the latitude/longitude coordinates corresponding to the center of the geoid and this will introduce an error in the conversion from latitude/longitude to geoid and back. For geoids at level 0, emprirical testing showed this to be less than 9.75 millimeters, which is sufficient for the type of applications intended for this libary. Note that this error is not constant accross the globe, and 9.75 millimeters is the maximum error seen.
procedure
(lat-lng-rect geoid) →
real? real? real? real? geoid : exact-integer?
The bounding box encloses the geoid minimally, but geoids and bounding boxes don’t overlap exactly, so the bounding box will always be bigger than then geoid and there will be other geoids which are inside the bounding box, but not in the geoid.
2 Storing and Retrieving from a SQLite database
procedure
(geoid->sqlite-integer geoid) → exact-integer?
geoid : exact-integer? (sqlite-integer->geoid i) → exact-integer? i : exact-integer?
Geoids are unsigned 64 bit values, and more than half of them have the highest bit set. SQLite will store numbers as signed 64 bits and convert unsigned 64 bit numbers greater than the maximum signed 64 bit value to floating point numbers, loosing precision. This means that geoids cannot be stored directly into a SQLite database.
These pair of functions subtract 263 from the geoid (or add that value back) to make sure the value is stored correctly. The ordering of the geoids is preserved, so they can still be used to index geograpic information.
3 Generating Test Data
procedure
(random-geoid level #:parent parent) → exact-integer?
level : exact-integer? parent : (or/c exact-integer #f)
This function is intended to generate geoids for unit tests.
procedure
(leaf-outline geoid [ #:steps steps #:closed? closed?]) → (listof exact-integer?) geoid : exact-integer? steps : (and/c integer? positive?) = 10 closed? : boolean? = #t
steps specifies the number of geoids ot put on each side of the rectangle, while closed? specifies if the first geoid should be duplicated as the last element in the list, to close the loop.
As with leaf-corners, geoids are placed in counter-clockwise order, but there is no guarantee about the start corner of the loop.
4 API Documentation
procedure
(last-valid-geoid) → exact-integer?
See geoid-stride to determine the amount to add to a geoid to obtain the next geoid at the same level.
procedure
This can be used to create half-open geoid ranges, for example by leaf-span.
procedure
(valid-geoid? geoid) → boolean?
geoid : exact-integer?
procedure
(geoid-level geoid) → (integer-in 0 30)
geoid : exact-integer?
Geoids produced by lat-lng->geoid are at level 0 and you can use enclosing-geoid to obtain a geoid at a higher level.
procedure
(geoid-stride geoid) → exact-integer?
geoid : exact-integer?
procedure
(enclosing-geoid geoid level) → exact-integer?
geoid : exact-integer? level : (integer-in 0 30)
procedure
(split-geoid geoid)
→ (list/c exact-integer? exact-integer? exact-integer? exact-integer?) geoid : exact-integer?
The returned geoids are not in any particular order.
procedure
(leaf-geoid? geoid) → boolean?
geoid : exact-integer?
procedure
(leaf-span geoid) →
exact-integer? exact-integer? geoid : exact-integer?
All geoids which are inside this geoid, regardless of level, are contained in the returned number range, so this range can be used to check if any geoid is inside this one.
The leaf span returned by this function can be used to search for geoids in an SQL query, however, if you do that, make sure you adjust them, as noted in the SQLite section above.
procedure
(leaf-span* geoids)
→ (list-of (cons/c exact-integer? exact-integer?)) geoids : (list-of exact-integer?)
This function can be used to create more efficient queries if a geoid is inside a list of geoids. A common use case is to use adjacent-geoids to obtain the neighbours of a geoid and using leaf-span* to find the ranges for all geoids in this neighbourhood.
procedure
(contains-geoid? this-geoid other-geoid) → boolean?
this-geoid : exact-integer? other-geoid : exact-integer?
This a convenient function, but if you need to check lots of geoids, this will be slower than obtainging the leaf-span of this-geoid and checking of other geoids are inside the returned range.
procedure
(leaf-corners geoid)
→ (list/c exact-integer? exact-integer? exact-integer? exact-integer?) geoid : exact-integer?
procedure
(adjacent-geoids geoid) → (list-of exact-integer?)
geoid : exact-integer?
Normally, 8 geoids are returned, but only 7 are returned if geoid is in a corner of a face and only 4 geoids if it is a face level geoid.
procedure
(approximate-area-for-geoid geoid) → real?
geoid : exact-integer? (approximate-area-for-level level) → real? level : (between/c 0 30)
These functions can be used to get a general idea of the surface covered by a geoid.
procedure
(distance-between-geoids g1 g2) → real?
g1 : exact-integer? g2 : exact-integer?
procedure
(distance-from-geoid g) → (-> exact-integer? real?)
g : exact-integer?
If you need to calculate the distance between a single geoid and several others, it might be faster to use this function to construct a "distance function".
5 Waypoint Alignment
(require geoid/waypoint-alignment) | package: geoid |
procedure
(waypoint-alignment-cost path1 path2) → real?
path1 : (vectorof exact-integer?) path2 : (vectorof exact-integer?)
A smaller cost indicates that the two paths are more similar. Ideally, the cost of a path against itself should be zero, but, due to floating point errors, this is a small positive number.
This function returns the Dynamic Time Warping cost of the two paths.
NOTE: the alignment cost will depend not only on how close the two paths are to each other, but also on the length of the paths, so it is up to the caller of the funtion to decide how to interpret the resulting cost and determine if the two paths are the same or not.
6 Region Tiling
(require geoid/tiling) | package: geoid |
The geoid/tiling module contains functions to calculate geoid coverings for regions covering the earth surface. A geoid covering is a list of geoids which are inside a region. This functionality can be used to implement fast region containment tests for geographic data.
procedure
(make-spherical-cap lat lon radius) → (is-a?/c region%)
lat : real? lon : real? radius : real?
This is not strictly a "region", but sometimes it is useful to determine the geoid covering for a stretch of road, and this can be done by defining the road as a sequence of points and tiling it with geoids.
procedure
(make-closed-polyline track #:ccw? ccw?) → (is-a?/c region%)
track : (listof (vector real? real?)) ccw? : #t
On a sphere (which the geoid library uses as a model for Earth), each closed list of points defines two regions, since the entire sphere surface is finite. The simplest way to visualize this is with a sequence of points along the Earths equator, which defines two regions: the northen hemisphere and the southern hemisphere. The ccw? parameter specifies which region is the "inside" region defined by the track: when ccw? is #t, the inside is the region which has the points in the counter-clockwise order, or to the left of the segments defined by the points as we go around the track.
See also guess-winding-order for determining the winding order of a sequence of points.
The closed polyline defined by track cannot have segments that intersect each other.
procedure
(join-regions r ...+) → (is-a?/c region%)
r : (is-a?/c region%)
procedure
(intersect-regions r ...+) → (is-a?/c region%)
r : (is-a?/c region%)
This function can be used to define regions from GeoJSON "polygon with holes" inputs – since the inner polygons have oposite winding order from the outer polygon, they can be created with the same winding order as the outer one, thus defining the opposite, or outer region, and intersecting them with the outer region.
Regions deifned in GeoJSON objects don’t specify a winding order, so this function can be used to determine the winding order of polygons inside GeoJSON objects, which can be passed on to make-closed-polyline.
procedure
(geoid-tiling-for-region r min-level max-level) → (listof integer?) r : (is-a?/c region%) min-level : integer? max-level : integer?
See geoid-level for a discution on the geoid levels.
This function calls coarse-geoid-tiling-for-region than refine-geoid-tiling-for-region for each geoid which intersects the region, and returns all the geoids in a single list.
procedure
(coarse-geoid-tiling-for-region r level)
→
(listof integer?) (listof integer?) r : (is-a?/c region%) level : integer?
See geoid-level for a discution on the geoid levels.
procedure
(refine-geoid-tiling-for-region r geoid min-level) → (listof integer?) r : (is-a?/c region%) geoid : integer? min-level : integer?
See geoid-level for a discution on the geoid levels.
7 Geodesy Calculations
(require geoid/geodesy) | package: geoid |
The geoid/geodesy module contains functions to calculate distances and bearings between points on the Earth surface. Two calculation modes are provided: one that approximates Earth as an ellipsoid, by default wgs84, which is the model used by the GPS satelites, and another which approximates the Earth as a sphere.
procedure
(ellipsoid? e) → boolean?
e : any/c
procedure
(make-ellipsoid major minor) → ellipsoid?
major : real? minor : real?
value
parameter
(geodesy-angle-mode) → (or/c 'degrees 'radians)
(geodesy-angle-mode mode) → void? mode : (or/c 'degrees 'radians)
= 'degrees
parameter
(geodesy-ellipsoid) → (or/c #f ellipsoid?)
(geodesy-ellipsoid ellipsoid) → void? ellipsoid : (or/c #f ellipsoid?)
= wgs84
procedure
(distance-between lat1 lon1 lat2 lon2 [ #:angle-mode m #:ellipsoid e]) → real? lat1 : real? lon1 : real? lat2 : real? lon2 : real? m : (or/c 'degrees 'radians) = (geodesy-angle-mode) e : (or/c #f ellipsoid?) = (geodesy-ellipsoid)
m specifies if the latitude and longitude angles are specified in degrees or radians, while e specifies the ellipsoid to use for the calculation. If e is #f, the calculation is done assuming Earth is a sphere.
procedure
(initial-bearing lat1 lon1 lat2 lon2 [ #:angle-mode m #:ellipsoid e]) → real? lat1 : real? lon1 : real? lat2 : real? lon2 : real? m : (or/c 'degrees 'radians) = (geodesy-angle-mode) e : (or/c #f ellipsoid?) = (geodesy-ellipsoid)
Note that, when traveling along a great circle the bearing will not remain constant.
The m and e parameters are the same as for distance-between.
procedure
(final-bearing lat1 lon1 lat2 lon2 [ #:angle-mode m #:ellipsoid e]) → real? lat1 : real? lon1 : real? lat2 : real? lon2 : real? m : (or/c 'degrees 'radians) = (geodesy-angle-mode) e : (or/c #f ellipsoid?) = (geodesy-ellipsoid)
procedure
(destination-point lat lon bearing distance [ #:angle-mode m #:ellipsoid e]) →
real? real? lat : real? lon : real? bearing : real? distance : real? m : (or/c 'degrees 'radians) = (geodesy-angle-mode) e : (or/c #f ellipsoid?) = (geodesy-ellipsoid)
The m and e parameters are the same as for distance-between.
procedure
(midway-point lat1 lon1 lat2 lon2 [ #:angle-mode m #:ellipsoid e]) →
real? real? lat1 : real? lon1 : real? lat2 : real? lon2 : real? m : (or/c 'degrees 'radians) = (geodesy-angle-mode) e : (or/c #f ellipsoid?) = (geodesy-ellipsoid)
The m and e parameters are the same as for distance-between.