1 /***********************************************************************
2  * Copyright (c) 2013-2024 Commonwealth Computer Research, Inc.
3  * All rights reserved. This program and the accompanying materials
4  * are made available under the terms of the Apache License, Version 2.0
5  * which accompanies this distribution and is available at
6  * http://www.opensource.org/licenses/apache2.0.php.
7  ***********************************************************************/
8 
9 package org.locationtech.geomesa.utils.geohash
10 
11 import com.typesafe.scalalogging.LazyLogging
12 import org.locationtech.jts.geom.{Coordinate, GeometryFactory, Point, PrecisionModel}
13 
14 import scala.collection.BitSet
15 import scala.collection.immutable.{BitSet => IBitSet}
16 
17 /**
18  * GeoHashes above GeoHash.MAX_PRECISION are not supported.
19  * @param x
20  * @param y
21  * @param bbox
22  * @param bitset
23  * @param prec
24  */
25 case class GeoHash private(x: Double,
26                            y: Double,
27                            bbox: BoundingBox,
28                            bitset: BitSet,
29                            prec: Int, // checked in factory methods in companion object
30                            private val optHash: Option[String]) extends Comparable[GeoHash] {
31   require(x >= -180.0 && x <= 180.0)
32   require(y >= -90.0  && y <= 90.0)
33 
34   import org.locationtech.geomesa.utils.geohash.GeoHash._
35 
36   /**
37    * Hash string is calculated lazily if GeoHash object was created
38    * from a Point, because calculation is expensive
39    */
40   lazy val hash = optHash.getOrElse(toBase32(bitset, prec))
41 
42   lazy val geom = bbox.geom
43 
44   /**
45    * Utility method to return the bit-string as a full binary string.
46    *
47    * NB:  We align our bits to the RIGHT by default, not to the LEFT.
48    * This means that the most significant bit in the GeoHash (as it is
49    * interpreted to a geometry) is actually the least significant bit
50    * in the bit-set.  In other words:  Interpret the bit-set from
51    * right to left.
52    *
53    * For clarity (?), this routine prints out the bits in MSB -> LSB
54    * order, so that their correspondence with the base-32
55    * characters is directly readable.
56    */
57   def toBinaryString: String =
58     (0 until prec).map(bitIndex => boolMap(bitset(bitIndex))).mkString
59 
60   /**
61     * Converts the geohash into an int.
62     *
63     * Note: geohashes with more than 6 base32 digits will not fit in an int
64     *
65     * @return the geohash as an int
66     */
67   def toInt: Int = {
68     var i = 0
69     var asInt = 0
70     while (i < prec) {
71       if (bitset(i)) {
72         asInt += 1 << prec - 1 - i
73       }
74       i += 1
75     }
76     asInt
77   }
78 
79   def getPoint = GeoHash.factory.createPoint(new Coordinate(x,y))
80 
81   def contains(gh: GeoHash): Boolean = prec <= gh.prec && bitset.subsetOf(gh.bitset)
82 
83   def next(): GeoHash =  GeoHash(GeoHash.next(bitset, prec), prec)
84 
85   override def equals(obj: Any): Boolean = obj match {
86     case that: GeoHash => this.bitset == that.bitset && this.prec == that.prec
87     case _ => false
88   }
89 
90   // Overriding equals obligates us to override hashCode.
91   override def hashCode: Int = bitset.hashCode + prec
92 
93   override def compareTo(gh: GeoHash) = this.hash.compareTo(gh.hash)
94 }
95 
96 case class Bounds(low: Double,
97                   high: Double) {
98   lazy val mid = (low+high)/2.0
99 }
100 
101 object GeoHash extends LazyLogging {
102 
103   val MAX_PRECISION = 63 // our bitset operations assume all bits fit in one Long
104 
105   implicit def toGeometry(gh: GeoHash) = gh.geom
106 
107   private val boolMap : Map[Boolean,String] = Map(false -> "0", true -> "1")
108 
109   lazy val factory: GeometryFactory = new GeometryFactory(new PrecisionModel, 4326)
110 
111   private val latBounds = Bounds(-90.0,90.0)
112   private lazy val latRange: Double = latBounds.high - latBounds.low
113   private val lonBounds = Bounds(-180.0,180.0)
114   private lazy val lonRange: Double = lonBounds.high - lonBounds.low
115 
116   private lazy val powersOf2Map: Map[Int, Long] =
117     (0 to   MAX_PRECISION).map(i => (i, 1L << i)).toMap // 1L << i == math.pow(2,i).toLong
118 
119   private lazy val latDeltaMap: Map[Int, Double]  =
120     (0 to MAX_PRECISION).map(i => (i, latRange / powersOf2Map(i / 2))).toMap
121   def latitudeDeltaForPrecision(prec: Int): Double = {
122     checkPrecision(prec)
123     latDeltaMap(prec)
124   }
125 
126   private lazy val lonDeltaMap: Map[Int, Double] =
127     (0 to MAX_PRECISION).map(i => (i, lonRange / powersOf2Map(i / 2 + i % 2))).toMap
128   def longitudeDeltaForPrecision(prec: Int): Double = {
129     checkPrecision(prec)
130     lonDeltaMap(prec)
131   }
132 
133   private val bits = Array(16,8,4,2,1)
134   protected[geohash] val base32 = "0123456789bcdefghjkmnpqrstuvwxyz"
135   private val characterMap: Map[Char, BitSet] =
136     base32.zipWithIndex.map { case (c, i) => c -> bitSetFromBase32Character(i) }.toMap
137 
138   private val lonMax = 360 - math.pow(0.5, 32)
139   private val latMax = 180 - math.pow(0.5, 32)
140 
141   // create a new GeoHash from a binary string in MSB -> LSB format
142   // (analogous to what is output by the "toBinaryString" method)
143   def fromBinaryString(bitsString: String): GeoHash = {
144     val numBits = bitsString.length
145     val bitSet: BitSet = bitsString.zipWithIndex.foldLeft(BitSet())((bsSoFar, bitPosTuple) => bitPosTuple match {
146       case (bitChar, pos) if bitChar == '1' => bsSoFar + pos
147       case _                                => bsSoFar
148     })
149     apply(bitSet, numBits)
150   }
151 
152   def apply(string: String): GeoHash = decode(string) // precision checked in decode
153   def apply(string: String, precision:Int): GeoHash = decode(string, Some[Int](precision)) // precision checked in decode
154 
155   // We expect points in x,y order, i.e., longitude first.
156   def apply(p: Point, prec: Int): GeoHash = apply(p.getX, p.getY, prec) // precision checked in apply
157   def apply(bs: BitSet, prec: Int): GeoHash = decode(toBase32(bs, prec), Some(prec)) // precision checked in decode
158 
159   // We expect points x,y i.e., lon-lat
160   def apply(lon: Double, lat: Double, prec: Int = 25): GeoHash = {
161     require(lon >= -180.0 && lon <= 180.0)
162     require(lat >= -90.0  && lat <= 90.0)
163     checkPrecision(prec)
164 
165     val lonDelta = lonDeltaMap(prec)
166     val latDelta = latDeltaMap(prec)
167     val lonIndex = if(lon == 180.0) (lonMax / lonDelta).toLong else ((lon - lonBounds.low) / lonDelta).toLong
168     val latIndex = if(lat == 90.0)  (latMax / latDelta).toLong else ((lat - latBounds.low) / latDelta).toLong
169 
170     encode(lonIndex, latIndex, lonDelta, latDelta, prec)
171   }
172 
173   def covering(ll: GeoHash, ur: GeoHash, prec: Int = 25) = {
174     checkPrecision(prec)
175 
176     val bbox = BoundingBox(ll.getPoint, ur.getPoint)
177 
178     def subsIntersecting(hash: GeoHash): Seq[GeoHash] = {
179       if (hash.prec == prec) List(hash)
180       else {
181         val subs = BoundingBox.generateSubGeoHashes(hash)
182         subs.flatMap { case sub =>
183           if (!bbox.intersects(sub.bbox)) List()
184           else subsIntersecting(sub)
185         }
186       }
187     }
188     val init = BoundingBox.getCoveringGeoHash(bbox, prec)
189     subsIntersecting(init)
190   }
191 
192   def checkPrecision(precision: Int) =
193     require(precision <= MAX_PRECISION,
194             s"GeoHash precision of $precision requested, but precisions above $MAX_PRECISION are not supported")
195 
196   /**
197    * Get the dimensions of the geohash grid bounded by ll and ur at precision.
198    * @param g1
199    * @param g2
200    * @param precision
201    * @return tuple containing (latitude span count, longitude span count)
202    */
203   def getLatitudeLongitudeSpanCount(g1: GeoHash, g2: GeoHash, precision: Int): (Int, Int) = {
204     require(g1.prec == precision,
205             s"Geohash ${g1.hash} has precision ${g1.prec} but precision ${precision} is required")
206     require(g2.prec == precision,
207             s"Geohash ${g2.hash} has precision ${g2.prec} but precision ${precision} is required")
208 
209     val Array(latIndex1, lonIndex1) = GeoHash.gridIndicesForLatLong(g1)
210     val Array(latIndex2, lonIndex2) = GeoHash.gridIndicesForLatLong(g2)
211 
212     ((math.abs(latIndex2 - latIndex1) + 1).toInt, (math.abs(lonIndex2 - lonIndex1) + 1).toInt)
213   }
214 
215   /**
216    * Convenience method to return both the latitude and longitude indices within
217    * the grid of geohashes at the precision of the specified geohash
218    *
219    * @param gh the geohash
220    * @return an array containing first the latitude index and then the longitude index
221    */
222   def gridIndicesForLatLong(gh: GeoHash) = {
223     val (lonIndex, latIndex) = extractReverseBits(gh.bitset.toBitMask.head, gh.prec)
224     Array(latIndex, lonIndex)
225   }
226 
227   /**
228    * Gets a long value representing a latitude index within the grid of geohashes
229    * at the precision of the specified geohash
230    *
231    * @param gh the geohash
232    * @return latitude index
233    */
234   def gridIndexForLatitude(gh: GeoHash) = {
235     val (_, latIndex) = extractReverseBits(gh.bitset.toBitMask.head, gh.prec)
236     latIndex
237   }
238 
239 
240   /**
241    * Gets a long value representing a longitude index within th grid of geohashes
242    * at the precision of the specified geohash
243    *
244    * @param gh the geohash
245    * @return longitude index
246    */
247   def gridIndexForLongitude(gh: GeoHash) = {
248     val (lonIndex, _) = extractReverseBits(gh.bitset.toBitMask.head, gh.prec)
249     lonIndex
250   }
251 
252 
253   /**
254    * Composes a geohash from a latitude and longitude index for the grid of geohashes
255    * at the specified precision.
256    *
257    * Note that the maximum latitude index is 2^(prec / 2) - 1 and the maximum longitude index
258    * is 2^(prec / 2 + prec % 2) -1 for the given precision.  An exception will be thrown if a
259    * larger index value is passed to this function.
260    *
261    * @param latIndex latitude index
262    * @param lonIndex longitude index
263    * @param prec the precision
264    * @return a geohash at the specified latitude and longitude index for the given geohash precision
265    */
266   def composeGeoHashFromBitIndicesAndPrec(latIndex: Long, lonIndex: Long, prec: Int): GeoHash = {
267     checkPrecision(prec)
268     encode(lonIndex, latIndex, lonDeltaMap(prec), latDeltaMap(prec), prec)
269   }
270 
271 
272   def next(gh: GeoHash): GeoHash = GeoHash(GeoHash.next(gh.bitset), gh.prec)
273 
274   def next(bs:BitSet, precision:Int=63) : BitSet = {
275     (0 until precision).reverse.foldLeft(true,BitSet())((t,idx) => t match { case (carry,newBS) => {
276       if (carry) {
277         if (bs(idx)) (true, newBS)
278         else (false, newBS+idx)
279       } else {
280         if (bs(idx)) (false, newBS+idx)
281         else (false, newBS)
282       }
283     }})._2
284   }
285 
286 
287   /**
288    * Create the geohash at the given latitude and longitude index with the given deltas and precision.
289    * Assumes prec <= 63, that is, all bits in the bitset fit in one Long
290    * @param lonIndex: the longitude index (x index)
291    * @param latIndex: the latitude index (y index)
292    * @param lonDelta the longitude (x) delta per grid cell for the given precision
293    * @param latDelta the latitude (y) delta per grid cell for the given precision
294    * @param prec precision (# of bits)
295    * @return the encoded GeoHash (note that calculation of the actual hash string is lazy)
296    */
297   def encode(lonIndex: Long, latIndex: Long, lonDelta: Double, latDelta: Double, prec: Int): GeoHash = {
298 
299     val bitSet = IBitSet.fromBitMaskNoCopy(Array(interleaveReverseBits(lonIndex, latIndex, prec)))
300     val bbox = bboxFromIndiciesDeltas(lonIndex, latIndex, lonDelta, latDelta)
301 
302     GeoHash(bbox.midLon, bbox.midLat, bbox, bitSet, prec, None)
303   }
304 
305   /**
306    * There is no visible difference between "t4bt" as a 20-bit GeoHash and
307    * the same string as a 17-bit GeoHash.  Unless otherwise specified, assume
308    * that the string represents a full complement of bits.
309    *
310    * If the call specifies a precision, then there is some additional work:
311    * 1.  the full-characters (those each representing a full 5 bits) should
312    *     be interpreted as they were before
313    * 2.  the partial-character (zero or one; representing the remainder bits,
314    *     guaranteed to be fewer than 5) must be appended
315    *
316    * @param string the base-32 encoded string to decode
317    * @param precisionOption the desired precision as an option
318    * @return the decoded GeoHash
319    */
320   private def decode(string: String, precisionOption: Option[Int] = None): GeoHash = {
321     // figure out what precision we should use
322     val prec = precisionOption.getOrElse(5*string.length)
323     checkPrecision(prec)
324 
325     // compute bit-sets for both the full and partial characters
326     val bitsets = string.zipWithIndex.map { case (c: Char, i: Int) => shift(i*5, characterMap(c)) }
327 
328     // OR all of these bit-sets together
329     val finalBitset: BitSet = bitsets.size match {
330       case 0 => BitSet()
331       case 1 => bitsets(0)
332       case _ => bitsets.reduce(_|_)
333     }
334 
335     val (lonIndex, latIndex) = extractReverseBits(finalBitset.toBitMask.head, prec)
336     val (lonDelta, latDelta) = (lonDeltaMap(prec), latDeltaMap(prec))
337     val bbox = bboxFromIndiciesDeltas(lonIndex, latIndex, lonDelta, latDelta)
338 
339     GeoHash(bbox.midLon, bbox.midLat, bbox, finalBitset, prec, Some(string))
340   }
341 
342   private def bboxFromIndiciesDeltas(lonIndex: Long,
343                                      latIndex: Long,
344                                      lonDelta: Double,
345                                      latDelta: Double): BoundingBox =
346     BoundingBox(Bounds((lonBounds.low+lonDelta*lonIndex), (lonBounds.low+lonDelta*(lonIndex+1))),
347                 Bounds((latBounds.low+latDelta*latIndex), (latBounds.low+latDelta*(latIndex+1))))
348 
349   /**
350    * Reverses and interleaves every other bits of two longs into one long. The two longs must be
351    * same size or first can be one bit longer than second.
352    * @param first first long
353    * @param second second long
354    * @param numBits The total number of bits of the interleaved & reversed result
355    * @return long with a total of numBits bits of first and second interleaved & reversed
356    */
357   private def interleaveReverseBits(first: Long, second: Long, numBits: Int): Long = {
358     /* We start with the first value of the interleaved long, coming from first if
359        numBits is odd or from second if numBits is even */
360     val even = (numBits & 0x01) == 0
361     val (actualFirst, actualSecond) = if(even) (second, first) else (first, second)
362     val numPairs = numBits >> 1
363     var result = 0L
364     (0 until numPairs).foreach{ pairNum =>
365       result = (result << 1) | ((actualFirst >> pairNum) & 1L)
366       result = (result << 1) | ((actualSecond >> pairNum) & 1L)
367     }
368     if (!even) (result << 1) | ((actualFirst >> numPairs) & 1L) else result
369   }
370 
371   /**
372    * Reverses and extracts every other bit up to numBits and extracts them into two longs.
373    * @param value the long
374    * @param numBits the number of bits to extract from the long
375    * @return (first long, second long)  The first long will be one bit longer than the second if
376    *         numBits is odd
377    */
378   private def extractReverseBits(value: Long, numBits: Int): (Long, Long) = {
379     val numPairs = numBits >> 1
380     var first = 0L
381     var second = 0L
382     (0 until numPairs).foreach{ pairNum =>
383       first = (first << 1) | ((value >> pairNum * 2) & 1L)
384       second = (second << 1) | ((value >> (pairNum * 2 + 1)) & 1L)
385     }
386     ((if((numBits & 0x01) == 1) ((first << 1) | ((value >> (numBits - 1)) & 1L)) else first), second)
387   }
388 
389   private def shift(n: Int, bs: BitSet): BitSet = bs.map(_ + n)
390 
391   private def bitSetFromBase32Character(charIndex: Long): BitSet =
392     BitSet(toPaddedBinaryString(charIndex, 5).zipWithIndex.filter { case (c,idx) => c == '1' }.map(_._2): _*)
393 
394   private def toPaddedBinaryString(i: Long, length: Int): String =
395     String.format("%" + length + "s", i.toBinaryString).replace(' ', '0')
396 
397   private def toBase32(bitset: BitSet, prec: Int): String = {
398     // compute the precision padded to the next 5-bit boundary
399     val numLeftoverBits = prec % 5
400     val precision = prec + (numLeftoverBits match {
401       case 0 => 0
402       case _ => (5-numLeftoverBits)
403     })
404 
405     // take the bit positions in groups of 5, and map each set to a character
406     // (based on the current bit-set); this works for off-5 precisions, because
407     // the additional bits will simply not be there (assumed to be zero)
408     (0 until precision).grouped(5).map(i=>ch(i, bitset)).mkString
409   }
410 
411   private def ch(v: IndexedSeq[Int], bitset: BitSet) =
412     base32(v.foldLeft(0)((cur,i) => cur + (if (bitset(i)) bits(i%bits.length) else 0)))
413 
414 }