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.geomesa.utils.geotools.GeometryUtils
13 import org.locationtech.geomesa.utils.iterators.CartesianProductIterable
14 import org.locationtech.geomesa.utils.text.WKTUtils
15 import org.locationtech.jts.geom._
16 import org.locationtech.spatial4j.context.jts.JtsSpatialContext
17 
18 import scala.collection.BitSet
19 import scala.collection.immutable.HashSet
20 import scala.collection.immutable.Range.Inclusive
21 import scala.util.Try
22 import scala.util.control.Exception.catching
23 
24 /**
25  * The following bits of code are related to common operations involving
26  * GeoHashes, such as recommending a GeoHash precision for an enclosing
27  * polygon; decomposing a polygon into a fixed number of subordinate
28  * GeoHashes; enumerating possible sub-strings within subordinate GeoHashes
29  * for a given polygon; etc.
30  */
31 object GeohashUtils
32   extends GeomDistance
33   with LazyLogging {
34 
35   // make sure the implicits related to distance are in-scope
36   import org.locationtech.geomesa.utils.geotools.Conversions.RichGeometry
37 
38   // these three constants are used in identifying unique GeoHash sub-strings
39   val Base32Padding = (0 to 7).map(i => List.fill(i)(GeoHash.base32.toSeq))
40   val BinaryPadding = (0 to 4).map(i => List.fill(i)(Seq('0', '1')))
41 
42   /**
43    * Simple place-holder for a pair of resolutions, minimum and maximum, along
44    * with an increment.
45    *
46    * @param minBitsResolution minimum number of bits (GeoHash) resolution to consider; must be ODD
47    * @param maxBitsResolution maximum number of bits (GeoHash) resolution to consider; must be ODD
48    */
49   case class ResolutionRange(minBitsResolution:Int=5,
50                              maxBitsResolution:Int=63,
51                              numBitsIncrement:Int=2) {
52     // validate resolution arguments
53     if (minBitsResolution >= maxBitsResolution)
54       throw new IllegalArgumentException("Minimum resolution must be strictly greater than maximum resolution.")
55 
56     lazy val range: Range = new Inclusive(minBitsResolution, maxBitsResolution, numBitsIncrement)
57     override def toString(): String = "{" + minBitsResolution.toString + ", +" + numBitsIncrement + "..., " + maxBitsResolution.toString + "}"
58 
59     def getNumChildren: Int = 1 << numBitsIncrement
60 
61     def getNextChildren(parent:BitSet, oldPrecision:Int) : List[BitSet] =
62       Range(0, getNumChildren).map { i =>
63         // compute bit-set corresponding to this integer,
64         // and add it to the left-shifted version of the parent
65         val bitString = i.toBinaryString
66 
67         Range(0, numBitsIncrement).foldLeft(parent) { case (bs, j) =>
68           val c = if (j < bitString.length) bitString.charAt(j) else '0'
69           if (c == '1') bs + (oldPrecision + bitString.length - 1 - j)
70           else bs
71         }
72       }.toList
73 
74   }
75 
76   // default precision model
77   val maxRealisticGeoHashPrecision : Int = 45
78   val numDistinctGridPoints: Long = 1L << ((maxRealisticGeoHashPrecision+1)/2).toLong
79   val defaultPrecisionModel = new PrecisionModel(numDistinctGridPoints.toDouble)
80 
81   // default factory for WGS84
82   val defaultGeometryFactory : GeometryFactory = new GeometryFactory(defaultPrecisionModel, 4326)
83 
84   /**
85    * Utility function that computes the minimum-bounding GeoHash that completely
86    * encloses the given (target) geometry.  This can be useful as a starting
87    * point for other calculations on the geometry.
88    *
89    * This is a bit of a short-cut, because what it *really* does is to start
90    * with the centroid of the target geometry, and step down through resolutions
91    * (in bits) for as long as a GeoHash constructed for that centroid continues
92    * to enclose the entire geometry.  This works, because any GeoHash that
93    * encloses the entire geometry must also include the centroid (as well as
94    * any other point).
95    *
96    * @param geom the target geometry that is to be enclosed
97    * @return the smallest GeoHash -- at an odd number of bits resolution -- that
98    *         completely encloses the target geometry
99    */
100   def getMinimumBoundingGeohash(geom:Geometry, resolutions:ResolutionRange) : GeoHash = {
101     // save yourself some effort by computing the geometry's centroid and envelope up front
102     val centroid = getCentroid(geom)
103     val env = defaultGeometryFactory.toGeometry(geom.getEnvelopeInternal)
104 
105     // conduct the search through the various candidate resolutions
106     val (_, ghOpt) = resolutions.range.foldRight((resolutions.minBitsResolution, Option.empty[GeoHash])){
107       case (bits, orig@(res, _)) =>
108         val gh = GeoHash(centroid.getX, centroid.getY, bits)
109         if (gh.contains(env) && bits >= res) (bits, Some(gh)) else orig
110     }
111 
112     // validate that you found a usable result
113     val gh = ghOpt.getOrElse(GeoHash(centroid.getX, centroid.getY, resolutions.minBitsResolution))
114     if (!gh.contains(env) && !gh.geom.equals(env))
115       throw new Exception("ERROR:  Could not find a suitable " +
116         resolutions.minBitsResolution + "-bit MBR for the target geometry:  " +
117         geom)
118 
119     gh
120   }
121 
122   /**
123    * Computes the centroid of the given geometry as a WGS84 point.
124    *
125    * @param geom the target geometry whose centroid is sought
126    * @return the centroid of the given geometry
127    */
128   def getCentroid(geom:Geometry) : Point = {
129     val pt = geom.safeCentroid()
130     geom.getFactory.createPoint(new Coordinate(pt.getX, pt.getY))
131   }
132 
133   // represents a degenerate (empty) geometry
134   lazy val emptyGeometry: Geometry = WKTUtils.read("POLYGON((0 0,0 0,0 0,0 0,0 0))")
135 
136   /**
137    * Utility class for the geometry-decomposition routines.  This represents
138    * one GeoHash-cell candidate within a decomposition.
139    *
140    * @param gh the associated GeoHash cell
141    * @param targetGeom the geometry being decomposed
142    */
143   abstract class DecompositionCandidate(val gh:GeoHash,
144                                         val targetGeom:Geometry,
145                                         val targetArea:Double,
146                                         val resolutions:ResolutionRange) {
147 
148     lazy val geomCatcher = catching(classOf[Exception])
149     lazy val area: Double = geomCatcher.opt{ gh.getArea }.getOrElse(0.0)
150     val areaOutside: Double
151     lazy val areaInside: Double = area - areaOutside
152     lazy val resolution: Int = gh.prec
153     lazy val isResolutionOk: Boolean =
154       resolution >= resolutions.minBitsResolution && resolution <= resolutions.maxBitsResolution
155 
156     lazy val intersectsTarget: Boolean =
157       geomCatcher.opt { gh.intersects(targetGeom) }.getOrElse(false)
158     lazy val intersection: Geometry =
159       geomCatcher.opt { gh.intersection(targetGeom) }.getOrElse(emptyGeometry)
160     lazy val intersectionArea: Double =
161       geomCatcher.opt { gh.intersection(targetGeom).getArea }.getOrElse(0.0)
162     def isLT(than:DecompositionCandidate): Boolean = {
163       if (areaOutside > than.areaOutside) true
164       else {
165         if (areaOutside == than.areaOutside) area < than.area
166         else false
167       }
168     }
169   }
170 
171   class PointDecompositionCandidate(gh:GeoHash,
172                                     targetGeom:Point,
173                                     targetArea:Double,
174                                     resolutions:ResolutionRange)
175     extends DecompositionCandidate(gh, targetGeom, targetArea, resolutions) {
176 
177     /**
178      * If the GeoHash does not contain the point, then the entire cell's area is
179      * outside of the target.  If the GeoHash does contain the point, then be
180      * careful:  Only some fraction of the cell's area should count as overage
181      * (otherwise, we can't favor smaller GeoHash cells in the decomposer).
182      */
183     override lazy val areaOutside: Double = area * (if (intersectsTarget) 0.75 else 1.0)
184   }
185 
186   class LineDecompositionCandidate(gh:GeoHash,
187                                    targetGeom:MultiLineString,
188                                    targetArea:Double,
189                                    resolutions:ResolutionRange)
190     extends DecompositionCandidate(gh, targetGeom, targetArea, resolutions) {
191 
192     /**
193      * If the GeoHash intersects the target lines, then the overlap is the
194      * area of the GeoHash cell less the length of the intersection.  Otherwise,
195      * they are disjoint, and the overlap is the entire area of the GeoHash cell.
196      *
197      * Yes, this mixes units, but it observes two trends:
198      * 1.  the longer a segment intersects, the smaller the area outside will be;
199      * 2.  the smaller a GeoHash cell, the smaller the area outside will be
200      */
201     override lazy val areaOutside : Double =
202       if (intersectsTarget) area * (1.0 - intersection.getLength / targetArea)
203       else area
204   }
205 
206   class PolygonDecompositionCandidate(gh:GeoHash,
207                                       targetGeom:MultiPolygon,
208                                       targetArea:Double,
209                                       resolutions:ResolutionRange)
210     extends DecompositionCandidate(gh, targetGeom, targetArea, resolutions) {
211 
212     /**
213      * If the GeoHash intersects the target polygon, then the overlap is the
214      * area of the GeoHash cell less the area of the intersection.  Otherwise,
215      * they are disjoint, and the overlap is the entire area of the GeoHash cell.
216      */
217     override lazy val areaOutside : Double =
218       if (intersectsTarget) area - intersection.getArea
219       else area
220   }
221 
222   def decompositionCandidateSorter(a:DecompositionCandidate,
223                                    b:DecompositionCandidate): Boolean = a.isLT(b)
224 
225   /**
226    * Decomposes the given polygon into a collection of disjoint GeoHash cells
227    * so that these constraints are satisfied:
228    * 1.  the total number of decomposed GeoHashes does not exceed the maximum
229    *     specified as an argument
230    * 2.  the resolution of the GeoHashes falls within the given range
231    * 3.  when replacing larger boxes with smaller boxes, always decompose first
232    *     the box that contains the most area outside the target polygon
233    *
234    * @param targetGeom the polygon to be decomposed
235    * @param maxSize the maximum number of GeoHash cells into which this polygon
236    *                should be decomposed
237    * @param resolutions the list of acceptable resolutions for the GeoHash cells
238    * @return the list of GeoHash cells into which this polygon was decomposed
239    *         under the given constraints
240    */
241   private def decomposeGeometry_(targetGeom: Geometry,
242                                  maxSize: Int = 100,
243                                  resolutions: ResolutionRange = new ResolutionRange(5,40,5)): List[GeoHash] = {
244     lazy val geomCatcher = catching(classOf[Exception])
245     val targetArea : Double = geomCatcher.opt { targetGeom.getArea }.getOrElse(0.0)
246     val targetLength : Double = geomCatcher.opt { targetGeom.getLength }.getOrElse(0.0)
247 
248     // qua factory
249     def createDecompositionCandidate(gh: GeoHash): DecompositionCandidate = {
250       // simple switch based on the geometry type
251       targetGeom match {
252         case multipoly: MultiPolygon    =>
253           new PolygonDecompositionCandidate(gh, multipoly, targetArea, resolutions)
254         case polygon: Polygon           =>
255           new PolygonDecompositionCandidate(
256             gh, new MultiPolygon(Array(polygon), polygon.getFactory), targetArea, resolutions)
257         case line: LineString           =>
258           new LineDecompositionCandidate(  // promote to a multi-line string of one element
259             gh, new MultiLineString(Array(line), line.getFactory), targetLength, resolutions)
260         case multiLine: MultiLineString =>
261           new LineDecompositionCandidate(gh, multiLine, targetLength, resolutions)
262         case point: Point               =>
263           new PointDecompositionCandidate(gh, point, targetArea, resolutions)  // should never be called, but it works
264         case _                          =>
265           throw new Exception(s"Unsupported Geometry type for decomposition:  ${targetGeom.getClass.getName}")
266       }
267     }
268 
269     // recursive routine that will do the actual decomposition
270     def decomposeStep(candidates: List[DecompositionCandidate]): List[DecompositionCandidate] = {
271       // complain, if needed
272       if (candidates.size > maxSize) throw new Exception("Too many candidates upon entry.")
273       else {
274         // identify the partial to replace...
275         // which of the single candidates contains the least overlap (area outside the target geometry)?
276         // assume that these are always sorted so that they are in descending order of overlap-area
277         val candidate : DecompositionCandidate = candidates(0)
278         val childResolution : Int = candidate.gh.prec+resolutions.numBitsIncrement
279 
280         // decompose this (worst) candidate into its four children
281         val candidateBitSet : BitSet = candidate.gh.bitset
282         val children:List[DecompositionCandidate] = resolutions.getNextChildren(candidateBitSet, candidate.gh.prec).map((childBitSet) => {
283           createDecompositionCandidate(GeoHash(childBitSet, childResolution))
284         }).filter(child => child.intersectsTarget)
285 
286         // build the next iteration of candidates
287         val newCandidates : List[DecompositionCandidate] = (candidates.tail ++ children).sortWith(decompositionCandidateSorter)
288 
289         // recurse, if appropriate
290         if ((newCandidates.size <= maxSize) && (childResolution <= resolutions.maxBitsResolution)) {
291           decomposeStep(newCandidates)
292         } else candidates
293       }
294     }
295 
296     // identify the smallest GeoHash that contains the entire polygon
297     val ghMBR = getMinimumBoundingGeohash(targetGeom, resolutions)
298     val candidateMBR = createDecompositionCandidate(ghMBR)
299 
300     // recursively decompose worst choices
301     val (keepers:List[DecompositionCandidate]) = decomposeStep(List(candidateMBR))
302 
303     // return only the keepers
304     (for (keeper:DecompositionCandidate <- keepers) yield keeper.gh).toList
305   }
306 
307   def getDecomposableGeometry(targetGeom: Geometry): Geometry = targetGeom match {
308     case g: Point                                            => targetGeom
309     case g: Polygon                                          => targetGeom
310     case g: LineString      if targetGeom.getNumPoints < 100 => targetGeom
311     case g: MultiLineString if targetGeom.getNumPoints < 100 => targetGeom
312     case _                                                   => targetGeom.convexHull
313   }
314 
315   /**
316    * Transforms a geometry with lon in (-inf, inf) and lat in [-90,90] to a geometry in whole earth BBOX
317    * 1) any coords of geometry outside lon [-180,180] are transformed to be within [-180,180]
318    *    (to avoid spatial4j validation errors)
319    * 2) use spatial4j to create a geometry with inferred International Date Line crossings
320    *    (if successive coordinates longitudinal difference is greater than 180)
321    * Parts of geometries with lat outside [-90,90] are ignored.
322    * To represent a geometry with successive coordinates having lon diff > 180 and not wrapping
323    * the IDL, you must insert a waypoint such that the difference is less than 180
324    */
325   def getInternationalDateLineSafeGeometry(targetGeom: Geometry): Try[Geometry] = {
326 
327     def degreesLonTranslation(lon: Double): Double = (((lon + 180) / 360.0).floor * -360).toInt
328 
329     def translateCoord(coord: Coordinate): Coordinate = {
330       new Coordinate(coord.x + degreesLonTranslation(coord.x), coord.y)
331     }
332 
333     def translatePolygon(geometry: Geometry): Geometry =
334       defaultGeometryFactory.createPolygon(geometry.getCoordinates.map(c => translateCoord(c)))
335 
336     def translateLineString(geometry: Geometry): Geometry =
337       defaultGeometryFactory.createLineString(geometry.getCoordinates.map(c => translateCoord(c)))
338 
339     def translateMultiLineString(geometry: Geometry): Geometry = {
340       val coords = (0 until geometry.getNumGeometries).map { i => geometry.getGeometryN(i) }
341       val translated = coords.map { c => translateLineString(c).asInstanceOf[LineString] }
342       defaultGeometryFactory.createMultiLineString(translated.toArray)
343     }
344 
345     def translateMultiPolygon(geometry: Geometry): Geometry = {
346       val coords = (0 until geometry.getNumGeometries).map { i => geometry.getGeometryN(i) }
347       val translated = coords.map { c => translatePolygon(c).asInstanceOf[Polygon] }
348       defaultGeometryFactory.createMultiPolygon(translated.toArray)
349     }
350 
351     def translateMultiPoint(geometry: Geometry): Geometry =
352       defaultGeometryFactory.createMultiPoint(geometry.getCoordinates.map(c => translateCoord(c)))
353 
354     def translatePoint(geometry: Geometry): Geometry = {
355       defaultGeometryFactory.createPoint(translateCoord(geometry.getCoordinate))
356     }
357 
358     def translateGeometry(geometry: Geometry): Geometry = {
359       geometry match {
360         case p: Polygon =>          translatePolygon(geometry)
361         case l: LineString =>       translateLineString(geometry)
362         case m: MultiLineString =>  translateMultiLineString(geometry)
363         case m: MultiPolygon =>     translateMultiPolygon(geometry)
364         case m: MultiPoint =>       translateMultiPoint(geometry)
365         case p: Point =>            translatePoint(geometry)
366       }
367     }
368 
369     Try {
370       // copy the geometry so that we don't modify the input - JTS mutates the geometry
371       // don't use the defaultGeometryFactory as it has limited precision
372       val copy = GeometryUtils.geoFactory.createGeometry(targetGeom)
373       val withinBoundsGeom =
374         if (targetGeom.getEnvelopeInternal.getMinX < -180 || targetGeom.getEnvelopeInternal.getMaxX > 180)
375           translateGeometry(copy)
376         else
377           copy
378 
379       JtsSpatialContext.GEO.makeShape(withinBoundsGeom, true, true).getGeom
380     }
381   }
382 
383   /**
384    * Quick-and-dirty sieve that ensures that we don't waste time decomposing
385    * single points.
386    */
387   def decomposeGeometry(targetGeom: Geometry,
388                         maxSize: Int = 100,
389                         resolutions: ResolutionRange = ResolutionRange(0, 40, 5),
390                         relaxFit: Boolean = true): List[GeoHash] = {
391   // quick hit to avoid wasting time for single points
392     targetGeom match {
393       case point: Point => List(GeoHash(point.getX, point.getY, resolutions.maxBitsResolution))
394       case gc: GeometryCollection => (0 until gc.getNumGeometries).toList.flatMap { i =>
395         decomposeGeometry(gc.getGeometryN(i), maxSize, resolutions, relaxFit)
396       }.distinct
397       case _ =>
398         val safeGeom = getInternationalDateLineSafeGeometry(targetGeom).getOrElse(targetGeom)
399         decomposeGeometry_(
400           if (relaxFit) getDecomposableGeometry(safeGeom)
401           else safeGeom, maxSize, resolutions)
402     }
403   }
404 
405   /**
406    * Given a collection of GeoHash hash sub-strings, this routine
407    * will build an iterator that generates all dotted variants.
408    * "Dotted" in this context means supporting the abstracted
409    * representation (in which higher-level hash-strings use periods
410    * for the base-32 characters they don't use).  For example,
411    * if the string is "dqb", then the dotted variants include all
412    * of the following:
413    *
414    *   dqb
415    *   dq.
416    *   d..
417    *   ...
418    *
419    * @param set the collection of GeoHash hash substrings to dot
420    * @param maxSize the maximum allowable number of entries in the final
421    *                iterator
422    * @return an iterator over the dotted collection of hash substrings,
423    *         constrained to one more than the maximum allowable size
424    *         (to enable overflow-detection on the outside)
425    */
426   def getGeohashStringDottingIterator(set: Seq[String], maxSize: Int): Iterator[String] = {
427     val len = set.headOption.map(_.length).getOrElse(0)
428     (for {
429       i <- (0 to len).iterator
430       hash <- set.map(_.take(i)).distinct
431       newStr = hash.take(i) + "".padTo(len - i, ".").mkString
432     } yield newStr).take(maxSize + 1)
433   }
434 
435   def promoteToRegion(geom: Geometry): Geometry = geom match {
436     case g: Point   =>
437       g.buffer(1e-6)
438     case g: Polygon =>
439       if (g.getArea > 0.0) g
440       else g.safeCentroid().buffer(1e-6)
441     case g          =>
442       val env = g.getEnvelope
443       if (env.getArea > 0.0) env
444       else env.getCentroid.buffer(1e-6)
445   }
446 
447   /**
448    * Given an index-schema format such as "%1,3#gh", it becomes necessary to
449    * identify which unique 3-character GeoHash sub-strings intersect the
450    * query polygon.  This routine performs exactly such an identification.
451    *
452    * The full GeoHashes from which the sub-strings are extracted are computed
453    * at 35 bits.
454    *
455    * Computing all of the 35-bit GeoHashes that intersect with the target
456    * geometry can take too long.  Instead, we start with the minimum-bounding
457    * GeoHash (which might be 0 bits), and recursively dividing it in two
458    * while remembering those GeoHashes that are completely contained in the
459    * target geometry.  This has a few advantages:
460    *
461    * 1.  we can stop recursing into GeoHashes at the coarsest
462    *     level (largest geometry) possible when they stop intersecting
463    *     the target geometry;
464    * 2.  instead of enumerating all of the GeoHashes that intersect, we
465    *     can stop as soon as we know that all possible children are known
466    *     to be inside the target geometry; that is, if a 13-bit GeoHash
467    *     is covered by the target, then we know that all 15-bit GeoHashes
468    *     that are its children will also be covered by the target
469    * 3.  if we ever find a GeoHash that is entirely covered by the target
470    *     geometry whose precision is no more than 5 times the "offset"
471    *     parameter's number of bits, then we can stop, because all possible
472    *     combinations are known to be used
473    *
474    * As an example, consider trying to enumerate the (3, 2) sub-strings of
475    * GeoHashes in a polygon that is only slightly inset within the
476    * Southern hemisphere.  This implicates a large number of 25-bit
477    * GeoHashes, but as soon as one of the GeoHashes that has 15 or fewer
478    * bits is found that is covered by the target, the search can stop
479    * for unique prefixes, because all of its 25-bit children will be
480    * distinct and will also be covered by the target.
481    *
482    * This is easier to explain with pictures.
483    *
484    * @param geom the query-polygon that must intersect candidate GeoHashes
485    * @param offset how many of the left-most GeoHash characters to skip
486    * @param length how many of the (remaining) GeoHash characters to use
487    * @param MAX_KEYS_IN_LIST the maximum allowable number of unique GeoHash
488    *                         sub-strings; when exceeded, the function returns
489    *                         an empty list
490    *
491    * @return the list of unique GeoHash sub-strings from 35-bits precision that
492    *         intersect the target polygon; an empty list if there are too many
493    */
494   def getUniqueGeohashSubstringsInPolygon(geom: Geometry,
495                                           offset: Int,
496                                           length: Int,
497                                           MAX_KEYS_IN_LIST: Int = Int.MaxValue - 1,
498                                           includeDots: Boolean = true): Try[Seq[String]] = Try {
499 
500     val cover = promoteToRegion(geom)
501 
502     //val cover: Geometry = geom.buffer(0)
503     val maxBits = (offset + length) * 5
504     val minBits = offset * 5
505     val usedBits = length * 5
506     val allResolutions = ResolutionRange(0, Math.min(35, maxBits), 1)
507     val maxKeys = Math.min(2 << Math.min(usedBits, 29), MAX_KEYS_IN_LIST)
508     val polyCentroid = cover.safeCentroid()
509 
510     // find the smallest GeoHash you can that covers the target geometry
511     val ghMBR = getMinimumBoundingGeohash(geom, allResolutions)
512 
513     // this case-class closes over properties of the current search
514     case class BitPrefixes(prefixes: Seq[String]) {
515 
516       val hasEverythingPrefix = prefixes.exists(prefix => prefix.length <= minBits)
517 
518       // how many GeoHashes are entailed by the list of prefixes
519       val entailedSize =
520         if (hasEverythingPrefix) maxKeys
521         else Math.min(
522           1 << usedBits,
523           prefixes.foldLeft(0)((sumSoFar, prefix) => {
524             sumSoFar + (1 << Math.min(usedBits, maxBits - prefix.length))
525           }))
526 
527       // is there any prefix wholly contained within the target geometry
528       // that uses fewer than 5*offset bits?  if so, then all possible
529       // sub-strings are entailed
530       val usesAll = prefixes.exists(prefix => prefix.length <= minBits) ||
531         entailedSize == maxKeys
532 
533       // the loose inequality is so that we can detect overflow
534       def hasRoom = entailedSize <= maxKeys
535 
536       def notDone = !usesAll && hasRoom
537 
538       def overflowed =
539         if (usesAll) {
540           (1 << usedBits) > maxKeys
541         } else {
542           entailedSize > maxKeys
543         }
544 
545       // generate all combinations of GeoHash strings of
546       // the desired length
547       def generateAll(prefix: String): Seq[String] = {
548         val prefixHash = GeoHash.fromBinaryString(prefix).hash
549         if (prefixHash.length < length) {
550           val charSeqs = Base32Padding(length - prefixHash.length)
551           CartesianProductIterable(charSeqs).toList.map(prefixHash + _.mkString)
552         } else Seq(prefixHash)
553       }
554 
555       // generate all of the combinations entailed by the prefixes identified,
556       // all of which will have bits that overlap with the requested substring
557       // (that is, if we want (3,2), then all of the prefixes we identified have
558       // somewhere between 15 and 25 bits in them)
559       //
560       // each prefix, then, needs to be expanded to a list of all combinations
561       // of bits that reach the next 5-bit boundary, and then those prefixes
562       // can be expanded to use all combinations of base-32 characters that
563       // allow them to fill out the requisite range
564       def generateSome: Seq[String] = {
565         prefixes.foldLeft(HashSet[String]())((ghsSoFar, prefix) => {
566           // fill out this prefix to the next 5-bit boundary
567           val bitsToBoundary = (65 - prefix.length) % 5
568           val bases =
569             if (bitsToBoundary == 0) Seq(prefix)
570             else {
571               val fillers = BinaryPadding(bitsToBoundary)
572               val result = CartesianProductIterable(fillers).toList.map(prefix + _.mkString)
573               result
574             }
575           bases.foldLeft(ghsSoFar)((ghs, base) => {
576             val baseTrimmed = base.drop(minBits)
577             val newSubs = generateAll(baseTrimmed)
578             ghs ++ newSubs
579           })
580         }).toSeq
581       }
582 
583       def toSeq: Seq[String] =
584         if (usesAll) generateAll("")
585         else generateSome
586     }
587 
588     // assume that this method is never called on a GeoHash
589     // whose binary-string encoding is too long
590     def considerCandidate(candidate: GeoHash): Seq[String] = {
591       val bitString = candidate.toBinaryString
592 
593       if (!geom.intersects(candidate.geom)) return Nil
594 
595       if (cover.covers(candidate.geom) || (bitString.size == maxBits)) {
596         Seq(bitString)
597       } else {
598         if (bitString.size < maxBits) {
599           // choose which direction to recurse into next by proximity
600           // of the two child GeoHashes to the polygon's centroid;
601           // for rectangles or polygons whose area is concentrated
602           // near the centroid, this provides for a a HUGE speed increase
603           val gh0 = GeoHash.fromBinaryString(bitString + "0")
604           val gh1 = GeoHash.fromBinaryString(bitString + "1")
605           val d0 = Math.hypot(gh0.getPoint.getX - polyCentroid.getX, gh0.getPoint.getY - polyCentroid.getY)
606           val d1 = Math.hypot(gh1.getPoint.getX - polyCentroid.getX, gh1.getPoint.getY - polyCentroid.getY)
607           val (firstChild, secondChild) =
608             if (d0 <= d1) (gh0, gh1)
609             else (gh1, gh0)
610 
611           val firstChildList = considerCandidate(firstChild)
612 
613           // if you've found an entry that entails all sub-strings, stop searching
614           firstChildList ++ (firstChildList.headOption match {
615             case Some(bitStr) if bitStr.length <= minBits => Nil
616             case _                                        =>
617               considerCandidate(secondChild)
618           })
619         } else Nil
620       }
621     }
622 
623     // compute the list of acceptable prefixes
624     val bitPrefixes = BitPrefixes(
625       if (ghMBR.prec <= maxBits) considerCandidate(ghMBR)
626       else Seq(ghMBR.toBinaryString.drop(minBits).take(usedBits)))
627 
628     // detect overflow
629     if (bitPrefixes.overflowed) throw new IllegalStateException("Bit prefixes overflowed while calculating unique Geohash substrings in polygon using the following parameters: " +
630       s"\nGeometry: $geom \nOffset: $offset \nLength: $length \nMax Keys in List: $MAX_KEYS_IN_LIST")
631 
632     // not having overflowed, turn the collection of disjoint prefixes
633     // into a list of full geohash substrings
634     val unDotted = bitPrefixes.toSeq
635 
636     // add dotted versions, if appropriate (to match decomposed GeoHashes that
637     // may be encoded at less than a full 35-bits precision)
638     if (includeDots) {
639       if (unDotted.size < maxKeys) {
640         // STOP as soon as you've exceeded the maximum allowable entries
641         val keepers = getGeohashStringDottingIterator(
642           unDotted, MAX_KEYS_IN_LIST).take(MAX_KEYS_IN_LIST + 1).toList
643         if (keepers.size <= MAX_KEYS_IN_LIST) keepers.toSeq else Seq()
644       } else Seq()
645     } else unDotted
646   }
647 }