2e:null 2f:["$","title",null,{"children":"Spherical point encoding | Jacob Rus | Observable"}] 34:T50f,md` Given a fixed budget of bits, how should we efficiently map points on a sphere (e.g. representing places on Earth) to binary numbers? The ideal would be for (a) every bit pattern to represent a distinct point on the sphere, with (b) those points in an equally spaced arrangement, (c) each point representing a uniform hexagonal area on the sphere, (d) calculated from conventional representations using a simple, computationally cheap, easy-to-reason-about formula, (e) organized in a regular grid, (f) represented in a consistent way across precision scales, so that adding extra zero bits leaves a point unchanged, and (g) leaving a resulting representation from which other desired calculations could be made simply. Unfortunately, these criteria are incompatible, and some of them are impossible to achieve. As a result, a wide variety of representations have been proposed and used over the years. The best representation to use depends on the use case. Some representations are good for raster images of the sphere, others for practical computation (e.g. reprojecting and plotting an outline map in real time or querying a spatial database), others for transmitting compressed geographical shape data, and traditional geographical coordinates are the most widespread and familiar. `35:Tc24,md` More recently, cardinal directions can be replaced with [signed][] numbers, and numbers use decimal (base ten) fractions, so [as long as we know the coordinate order][lonlat] we can rewrite that location as (−122.426856, 37.758736). [Sometimes][geobuf] such values are multiplied by one million and stored as binary (base two) integers. This yields 360 × 10⁶ possible values or ${Math.log2(10**6 * 360).toFixed(1)} bits for longitude, and 1 bit fewer for latitude. Since the earth has a circumference of about 40,000 kilometers, this lets us represent points at the equator to within about 11 centimeters in each dimension. [signed]: https://en.wikipedia.org/wiki/Sign_(mathematics) [lonlat]: https://macwright.org/lonlat/ [geobuf]: https://github.com/mapbox/geobuf/blob/master/geobuf.proto In modern computers, fast arithmetic is done with binary [floating point][float] numbers, but this representation is not terribly well matched to representing angle measures. (Floating point numbers are roughly [log][] spaced, so many pairs do not represent valid points on the sphere, and those that do provide non-uniform precision.) As an alternative we can use the [fixed-point][] [binary angle][] representation, which would be the most uniform/efficient representation mentioned so far. If we used 32-bit integers, we could get ${(32 - Math.log2(10**6 * 360)).toFixed(1)} bits more precision for longitude and ${(32 - Math.log2(10**6 * 180)).toFixed(1)} bits more precision for latitude, compared to millionths of a degree. However, it typically makes sense to use the same scale for both latitude and longitude at the equator; in that case, if the two fixed-point numbers use the same precision, we don’t need one of the bits in the latitude number, and can either leave it unused or interpret it in some other way. Since 2³² ≈ 4.3 billion, with 32 bits of longitude and 31 bits of latitude we can represent points at the equator to roughly the nearest centimeter in each dimension. [float]: https://en.wikipedia.org/wiki/IEEE_754 [log]: https://en.wikipedia.org/wiki/Logarithmic_scale [binary angle]: https://en.wikipedia.org/wiki/Binary_scaling#Binary_angles [fixed-point]: https://en.wikipedia.org/wiki/Fixed-point_arithmetic At high latitudes (i.e. near the poles) the circles of constant latitude shrink so the longitudes become more closely spaced. In the limit, as latitude approaches a quarter turn, the circle of constant latitude shrinks down to a point. Thus if we want to represent a complicated shape near the poles, we end up using more precision for longitude than for latitude. For most purposes this extra precision is wasted, so we don’t lose anything if we more aggressively round longitudes near the poles. If ${$`\Delta \lambda`} represents the least representable change of longitude, we want ${$`\Delta \lambda \cos{\phi}`} to remain within 1 bit of the precision at the equator as latitude changes. If we have a compression algorithm which can detect bit-level redundancy, this can even save bandwidth. `36:T4bb,// take in longitude and latitude as floats in the ranges [-180, 180) and [-90, 90] and // emit floats in the ranges [-1, 1 - 2^(-32)) and [-0.5, 0.5] // with 1 representing π radians = 180 degrees. // consistently round latitude to some desired precision, and also round longitude // as aggressively as possible while maintaining the precision of longitude quantize_lonlat = { const ROUND_TO_INT = 2**53; // to round to nearest integer return function (longitude, latitude, precision=2**-31) { const slat = Math.sign(latitude), slon = Math.sign(longitude); latitude /= 180 * slat; longitude /= 180 * slon; if (latitude > 0.5) latitude = 0.5; // clamp latitude const latcos = Math.cos(Math.PI * latitude), latround = prev_pow2(ROUND_TO_INT * precision), lonround = prev_pow2(ROUND_TO_INT * precision / latcos); latitude = (latitude - latround) + latround; // round latitude longitude = (longitude - lonround) + lonround; // round longitude latitude *= slat; // restore signs longitude *= slon; longitude -= 2 * Math.round(longitude * 0.5); // map longitude to standard interval return [longitude, latitude]; } }37:T790,Cartesian = md` ### Cartesian coordinates Using angle measures requires computationally expensive and fiddly [circular functions][] rather than basic arithmetic to do pretty much anything, e.g. calculating distances and areas, rotating points on the sphere, judging orientation of shapes, and so on. Basic operations become much simpler if we switch to a 3-dimensional [Cartesian coordinate system][], and indeed many operations on geographical coordinates more or less amount to converting to Cartesian coordinates, performing some operation, then converting back, sometimes with a few tricky shortcuts to save work. We could instead directly store Cartesian coordinates and save the round-trips. If we assume the earth to be a unit-diameter sphere, [circular functions]: https://en.wikipedia.org/wiki/Circular_functions [Cartesian coordinate system]: https://en.wikipedia.org/wiki/Cartesian_coordinate_system ${$$` \begin{pmatrix} \lambda \\ \phi \end{pmatrix} \mapsto \begin{pmatrix} \cos{\lambda}\cos{\phi} \\[3mu] \sin{\lambda}\cos{\phi} \\[3mu] \sin{\phi} \end{pmatrix} `} Then the point at the playground could be represented as (−0.423935478, −0.66732383, 0.61233783). Unfortunately, for a given bit budget, this wrecks our precision, because we now need 3 coordinates instead of two. One possible workaround would be to store the signs of the three coordinates, along with two flag bits indicating which coordinate has the largest magnitude, and then store magnitudes of the other two coordinates using fixed point numbers between ${$`0..2^{-1/2}`}. The largest magnitude can then be inferred because ${$`x^2 + y^2 + z^2 = 1`}. But this is inelegant and still not especially effective. I think it best to leave Cartesian coordinates as a possible high-precision internal representation for doing arithmetic, but set it aside when storing or transmitting large amounts of 2-dimensional point data. ` 38:T6b5,Stereographic = md` ### Stereographic coordinates Another way to get from three coordinates down to two, which matches pretty well with the floating point arithmetic used by computers is the [stereographic projection][]. The stereographic projection is a perspective projection through one pole of the sphere onto a perpendicular plane. It is conformal (preserves local angles) and maps circles on the sphere to circles on the plane. [stereographic projection]: https://en.wikipedia.org/wiki/Stereographic_projection ${$$` \begin{pmatrix} \lambda \\ \phi \end{pmatrix} \mapsto \begin{pmatrix} \tan{\tfrac{1}{2}\big(\tfrac{\vphantom{1}\pi}{2} - \phi\big)} \cos{\lambda} \\[7mu] \tan{\tfrac{1}{2}\big(\tfrac{\vphantom{1}\pi}{2} - \phi\big)} \sin{\lambda} \end{pmatrix} `} This representation has the advantage over geographical coordinates that it does not need to cut the globe along a meridian, but only puncture it at a single point. (And since floating point numbers can represent ∞, this is no problem.) Also, it makes many kinds of later spherical geometry computations much simpler than they would be using the transcendental [circular functions][] required for geographical coordinates. [circular functions]: https://en.wikipedia.org/wiki/Trigonometric_functions If we use 32-bit [projective numbers][pj] to represent our stereographic coordinates of the Earth, at the worst we have a precision of about 3 centimeters in one coordinate. [pj]: /@jrus/pj To avoid unnecessary precision, we should round both coordinates to keep the values of ${$`\Delta x / (1 + x^2 + y^2)`} and ${$`\Delta y / (1 + x^2 + y^2)`} to within about 1 bit across the sphere. `39:T45b,// take in stereographic coordinates as floats and emit floats which // have been rounded to consistently match some desired precision quantize_stereo = { const ROUND_TO_INT = 2**53; // to round to nearest integer return function (x, y, precision=2**-30) { const sx = Math.sign(x), sy = Math.sign(y); x *= sx, y *= sy; const w = 0.5 * (1 + x*x + y*y); if (w > 2/(precision*precision)) return [Infinity, Infinity]; const round = prev_pow2(ROUND_TO_INT * precision * w), foo = Math.log2(precision * w), x2 = (x - round) + round, y2 = (y - round) + round, // potentially round the original number again using the precision // from the latitude of the first rounding. this is probably not // necessary in practice, but could avoid a case where two points // just across a rounding threshhold cross each-other. w2 = Math.min(w, 0.5 * (1 + x2*x2 + y2*y2)), round2 = prev_pow2(ROUND_TO_INT * precision * w2); x = (x - round2) + round2; y = (y - round2) + round2; x *= sx, y *= sy; // restore signs return [x, y]; } }3a:T778,md` ### Coordinates based on polyhedra We can choose a polyhedron (e.g. a [cube][], [octahedron][], [icosahedron][], [rhombic triacontahedron][], or *n*-gonal [bipyramid](https://en.wikipedia.org/wiki/bipyramid)) and a [spherical polyhedron][] of the same structure and associate each polygonal face with a spherical polygon, come up with a way of [projecting][] one to the other (e.g. using the [gnomonic projection][] or a [conformal projection][], and then come up with a grid on each polygon. [cube]: https://en.wikipedia.org/wiki/cube [octahedron]: https://en.wikipedia.org/wiki/octahedron [icosahedron]: https://en.wikipedia.org/wiki/icosahedron [rhombic triacontahedron]: https://en.wikipedia.org/wiki/Rhombic_triacontahedron [spherical polyhedron]: https://en.wikipedia.org/wiki/Spherical_polyhedron [projecting]: https://en.wikipedia.org/wiki/Map_projection [gnomonic projection]: https://en.wikipedia.org/wiki/Gnomonic_projection [conformal projection]: https://en.wikipedia.org/wiki/Schwarz%E2%80%93Christoffel_mapping One of the simplest and computationally cheapest possibilities is to use a gnomonic projection from each octant of the globe onto each face of an octahedron, but then use hexagonal pixels on each equilateral-triangle face. If we start with a triangular grid with 2^*n* divisions passing through the vertices of the octahedron, we end up with awkward decisions to make about which octant each point on a boundary fits into. We can instead offset our grid by ⅓ away from the edges in 4 non-adjacent octants and ⅔ away from the edges in the other 4 octants. This gives us overall a grid with 2^*n*+1 by 2^*n*+1 vertices, which can be conveniently stored in a square texture, in a [quincunxial][] arrangement. [quincunxial]: https://en.wikipedia.org/wiki/Peirce_quincuncial_projection To convert from cartesian coordinates to `3b:[] 32:{"id":"5d7102892473ca50","slug":"spherepoints","title":"Spherical point encoding","thumbnail":"c1933d2335905f5b2311a542b2198c7952fa3dc9be093201286332e0d40e8699","access_level":"unlisted","create_time":"$D2019-11-08T22:17:47.834Z","update_time":"$D2021-03-20T05:33:40.755Z","view_time":null,"first_public_time":null,"has_unpublished_edits":null,"fork":false,"stars":0,"starred":false,"owner":{"type":"team","id":"9b8a8e177a49059a","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":"starter_2024"},"version":711,"publish_version":711,"role":null,"authored":false,"authors":[{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null}],"members":[{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":0,"collection":null,"type":"notebook","body":{"title":"Spherical point encoding","cells":[{"id":0,"value":"md`# Spherical point encoding\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":5,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":263,"value":"md`\n### Geographical coordinates\n\nThe traditional method, extracted from ancient Mesopotamian astronomy, is to use the angle measures longitude ${$`\\lambda`} and latitude ${$`\\phi`} (the so-called [geographic coordinate system][], represented as [sexagesimal][] (base sixty) fractions of a 360° circumference, measuring East–West from the [prime meridian][] and North–South from the [equator][].\nA precise spot in my local playground might be written as (122°25′36″41‴ W, 37°45′31″27‴ N).\n\n [geographic coordinate system]: https://en.wikipedia.org/wiki/Geographic_coordinate_system\n [sexagesimal]: https://en.wikipedia.org/wiki/Sexagesimal\n [prime meridian]: https://en.wikipedia.org/wiki/Prime_meridian\n [equator]: https://en.wikipedia.org/wiki/Equator\n\n`\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":290,"value":"{\n var λ = -122.426856;\n var ϕ = 37.758736;\n var z = coord_zoom;\n const map = osmMap(λ, ϕ, z, Math.min(width, 640), 250);\n map.style.setProperty('border', '1px solid black')\n return html`

${map}`;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":324,"value":"viewof coord_zoom = Range([5, 19], {value: 12, step: 1, label: \"Zoom\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":72,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":426,"value":"prev_pow2 = (x) =>\n 2**(Math.floor(Math.log2(x)))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":148,"value":"$36","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":67,"value":"$37","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":80,"value":"$38","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":503,"value":"stereo = (lon, lat) => {\n var r = Math.tan(Math.PI * (90 - lat) / 360);\n lon *= Math.PI / 180;\n return [r * Math.cos(lon), r * Math.sin(lon)];\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":533,"value":"a = quantize_stereo(...stereo(2, -89.9999), 2**-prec)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":553,"value":"binaryInput({\n type: 'float64',\n value: a[0],\n displayParts: true\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":561,"value":"binaryInput({\n type: 'float64',\n value: a[1],\n displayParts: true\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":564,"value":"viewof prec = Range([0, 31], {value: 20, label: 'Precision (bits)'})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":602,"value":"c = a.map(pj_with_nbits(32))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":608,"value":"binaryInput({\n type: 'int32',\n value: c[0]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":613,"value":"binaryInput({\n type: 'int32',\n value: c[1]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":507,"value":"$39","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":555,"value":"import {binaryInput} from '@rreusser/binary-input'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":604,"value":"import {pj_with_nbits} from '@jrus/pj'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":112,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":226,"value":"md`* * *`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":399,"value":"viewof rnd = Range([0, 30], {step:1})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":395,"value":"b = 53 + rnd","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":349,"value":"'0b0' + (((0x2aaaaaaa - 2**b ) + 2**b) | ((0x55555555 - 2**b ) + 2**b)).toString(2)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":224,"value":"world = fetch(\"https://cdn.jsdelivr.net/npm/world-atlas@2/land-50m.json\").then(response => response.json())","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":582,"value":"import {math_css, $, $$} with {$css as $css} from \"@jrus/misc\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":584,"value":"$$css = html`${math_css}`","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":266,"value":"import {osmMap} from '@mourner/simple-web-map@963'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":327,"value":"import {Range} from '@observablehq/inputs'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W3b","stdlib":"1","titling":"auto"}} 33:"$32:body"