geoid – Efficient storage and queriying of geographic data
1 Converting to and from Latitude/  Longitude
lat-lng->geoid
geoid->lat-lng
lat-lng-rect
2 Storing and Retrieving from a SQLite database
geoid->sqlite-integer
sqlite-integer->geoid
3 Generating Test Data
random-geoid
leaf-outline
4 API Documentation
first-valid-geoid
last-valid-geoid
sentinel-geoid
valid-geoid?
geoid-level
geoid-stride
enclosing-geoid
split-geoid
leaf-geoid?
leaf-span
leaf-span*
contains-geoid?
leaf-corners
adjacent-geoids
approximate-area-for-geoid
approximate-area-for-level
distance-between-geoids
distance-from-geoid
5 Waypoint Alignment
waypoint-alignment-cost
6 Region Tiling
region%
make-spherical-cap
make-open-polyline
make-closed-polyline
join-regions
intersect-regions
subtract-regions
guess-winding-order
geoid-tiling-for-region
coarse-geoid-tiling-for-region
refine-geoid-tiling-for-region
7 Geodesy Calculations
ellipsoid?
make-ellipsoid
wgs84
geodesy-angle-mode
geodesy-ellipsoid
distance-between
initial-bearing
final-bearing
destination-point
midway-point
8.12

geoid – Efficient storage and queriying of geographic data🔗ℹ

Alex Harsányi

 (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.

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?
Convert a geographic coordinate (latidude, longitude) to and from a geoid. lat-lng->geoid will return a geoid at the specified level of precision, by default at the highest level of precision, level 0. geoid->lat-lng will accept geoids at any level of precision.

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?
Return the latitude / longitude rectangle which encloses the geoid. The four values returned are: minimum latitude, minimum longitude, maximum latitude, maximum longitude. The bounds are slightly extended, to ensure all leaf geoids are inside the bounding box.

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🔗ℹ

Convert a geoid into an integer suitable for storage into a SQLite database, or convert an integer stored in the database back into a geoid.

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)
Generate a random geoid at the specified level. If parent is specified, the level has to be lower (mode detailed) than the parent and a geoid which is inside the parent will be generated.

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
Return a sequence of leaf geoids representing the outline of geoid. This can be used to obtain the outline of a geoid to display on a map.

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🔗ℹ

Return the first and last valid integers that represent geoids. All integers between these two numbers are valid geoids, but they can be at different levels. Note that 0 is not a valid geoid.

See geoid-stride to determine the amount to add to a geoid to obtain the next geoid at the same level.

Returns an integer which is one bigger than the largest valid geoid (as returned by last-valid-geoid). The returned value is not a valid geoid.

This can be used to create half-open geoid ranges, for example by leaf-span.

procedure

(valid-geoid? geoid)  boolean?

  geoid : exact-integer?
Return true if geoid is a valid geoid. This is the same as testing that the geoid is between first-valid-geoid and last-valid-geoid, but shorter to type.

procedure

(geoid-level geoid)  (integer-in 0 30)

  geoid : exact-integer?
Return the level of this geoid – geoids can represent areas at increasing level of detail. The highest resolution is at level 0, and going up, geoids become 4 times bigger at each level. Level 30 is the top level, where the entire Earth surface is divided into 6 faces.

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?
Return the integer step that must be added to a geoid to obtain another geoid at the same level. This can be used to generate valid geoids in sequence. Note that the geoids at the highest resolution (level 2) have a stride of 2, so you cannot simply incremenet geoids.

procedure

(enclosing-geoid geoid level)  exact-integer?

  geoid : exact-integer?
  level : (integer-in 0 30)
Return the geoid at level that contains the specified geoid. This function can be used to obtain the geoid at a lower resolution which contains a given point or geoid.

Return the four geoids at the lower level into which the current geoid can be split. An error is raised, if the supplied geoid is a leaf geoid, at level 0.

The returned geoids are not in any particular order.

procedure

(leaf-geoid? geoid)  boolean?

  geoid : exact-integer?
Return #t if geoid is a geoid at level 0 (highest level of detail). This is equivalent to checking if geoid-level returns 0, but faster.

procedure

(leaf-span geoid)  
exact-integer? exact-integer?
  geoid : exact-integer?
Returns the half-open geoid range that are valid geoids contained in geoid. The first value is the smallest leaf geoid which is inside this geoid, the second value is either the smallest leaf geoid which is not part of geoid, or the sentinel-geoid, whichever is smaller.

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?)
Return a list of half open geoid ranges which define the valid geoids contained in all the geoids from geoids list. This is equivalent to calling leaf-span for each individual geoid in geoids and merging all the adjacent ranges together.

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?
Return true if the other-geoid is geographically inside this-geoid.

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.

Return the four leaf geoids which represent the corners of this geoid. The corners are returned in couter-clockwise order, but there is no guarante of which one is first (i.e. there is no guarantee that the list will start with the top-left corner)

procedure

(adjacent-geoids geoid)  (list-of exact-integer?)

  geoid : exact-integer?
Return the adjacent geoids which border geoid. The returned geoids will be at the same level as geoids.

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)
Return the approximate area, in square meters, of geoid or level. The area is calculated by dividing the earth surface by the number of geoids at that level and it is approximate because there will be geoids with smaller area and larger area than this at each level.

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?
Return the distance, in meters, on the Earth surface between the center of geoids g1 and g2. For leaf geoids this is a good appoximation for the distance between the locations represented by these geoids, since the size of a leaf geoid is approximately 8.5 millimeters.

procedure

(distance-from-geoid g)  (-> exact-integer? real?)

  g : exact-integer?
Return a function which can be used to calculate the distance on the Earth surface between the geoid g and another geoid. The returned function will accept a single geoid as an argument and will return the distance between that geoid and g.

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?)
Return a number representing how similar path1 is to path2. Both paths are vector of integers representing geoids.

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.

class

region% : class?

  superclass: object%

Region objects define a region on the earth surface and are used for tiling, but otherwise they are opaque objects.

procedure

(make-spherical-cap lat lon radius)  (is-a?/c region%)

  lat : real?
  lon : real?
  radius : real?
Create a region% which is a circle around the lat/lon point, where radius specifies the radius of the circle, in meters. The circle extends around the sphere representing Earth, thus the name "spherical cap".

procedure

(make-open-polyline track)  (is-a?/c region%)

  track : (listof (vector real? real?))
Create a region% representing the lines connecting lat/lon points in track.

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
Create a region% representing the inside region defined by the lines connecting lat/lon points in track.

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%)
Create a new region which is the union of the regions passed in as input.

procedure

(intersect-regions r ...+)  (is-a?/c region%)

  r : (is-a?/c region%)
Create a new region which is the intersection of the regions passed in as input.

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.

procedure

(subtract-regions r1 r2)  (is-a?/c region%)

  r1 : (is-a?/c region%)
  r2 : (is-a?/c region%)
Create a region containg the points which are inside r1 but not inside r2.

procedure

(guess-winding-order track)  (or/c 'cw 'ccw)

  track : (listof (vector real? real?))
Guess the winding order (clockwise or counter-clockwise) for the polygon defined by track, which is a list of latitude/longitude points. The guess works for "usual" polygons which are much smaller than a hemisphere and assumes that the "inside" is the smaller part.

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?
Return a list of geoids which cover the region r. The region is covered with geoids of level max-level first, than refined along the region borders with geoids down to min-level.

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?
Return two lists of geoids at level which cover the region r. All geoids will be the same level. The first list are all the geoids which are entirely inside the region r, while the second list contains all the geoids which intersect the region.

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?
Return a list of geoids which cover the region r up to min-level obtained by recursively splitting the input geoid.

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
Return #t if e is an ellipsoid created by make-ellipsoid.

procedure

(make-ellipsoid major minor)  ellipsoid?

  major : real?
  minor : real?
Create an ellipsoid with the major and minor semi-axes. Ellipsoids can be used as values for geodesy-ellipsoid or to individual functions in this module. The wgs84 ellipsoid, used by the GPS satelites, is already defined in this module and it is the default.

value

wgs84 : ellipsoid?

The WGS84 ellipsoid.

parameter

(geodesy-angle-mode)  (or/c 'degrees 'radians)

(geodesy-angle-mode mode)  void?
  mode : (or/c 'degrees 'radians)
 = 'degrees
A parameter that specifies the type of angles which are passed to the functions in this module as well as the type of angle which are returned for bearings. Valid values are 'degees, which means that latitude, longitudes are in degrees or 'radians.

parameter

(geodesy-ellipsoid)  (or/c #f ellipsoid?)

(geodesy-ellipsoid ellipsoid)  void?
  ellipsoid : (or/c #f ellipsoid?)
 = wgs84
A parameter that specifies the ellipsoid to use for the calculations by the functions in this module. By default, the wgs84 ellipsoid is used and a value of #f means that the calculations are done assuming a sperical Earth model.

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)
Return the distance between two points on Earth along the great circle, which is the shortest distance. The points are identified by latitude and longitude, which are angles.

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)
Return the initial bearing for traveling between two points on Earth along the great circle. The points are identified by latitude and longitude, which are angles and the bearning is also an angle, where 0 is the North direction and the angle moves clockwise.

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)
Similar to initial-bearing, but return the final bearing for traveling between two points on Earth along the great circle.

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)
Return the destination point for traveling from the point lat, lon along a great circle with an initial bearing for a given distance.

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)
Return the latitude/longitude point that it half-way between the points lat1, lon1 and lat2, lon2, along the great circle arc, that is along the shorted distance path between the two points.

The m and e parameters are the same as for distance-between.