org.locationtech.spatial4j.distance.DistanceUtils Java Examples

/**
 * Convert degrees to another unit
 *
 * @param degs  degrees
 * @param units UoM IRI of the unit to convert to
 * @return converted value as a double
 * @throws ValueExprEvaluationException
 */
public static double convertFromDegrees(double degs, IRI units) throws ValueExprEvaluationException {
	double v;

	if (GEOF.UOM_DEGREE.equals(units)) {
		v = degs;
	} else if (GEOF.UOM_RADIAN.equals(units)) {
		v = DistanceUtils.toRadians(degs);
	} else if (GEOF.UOM_UNITY.equals(units)) {
		v = degs / 180.0;
	} else if (GEOF.UOM_METRE.equals(units)) {
		v = DistanceUtils.degrees2Dist(degs, DistanceUtils.EARTH_MEAN_RADIUS_KM) * 1000.0;
	} else {
		throw new ValueExprEvaluationException("Invalid unit of measurement: " + units);
	}
	return v;
}
 

private void _checkHits(boolean bbox, Point pt, double distKM, int assertNumFound, int... assertIds) {
  SpatialOperation op = SpatialOperation.Intersects;
  double distDEG = DistanceUtils.dist2Degrees(distKM, DistanceUtils.EARTH_MEAN_RADIUS_KM);
  Shape shape = shapeFactory.circle(pt, distDEG);
  if (bbox)
    shape = shape.getBoundingBox();

  SpatialArgs args = new SpatialArgs(op,shape);
  //args.setDistPrecision(0.025);
  Query query = strategy.makeQuery(args);
  SearchResults results = executeQuery(query, 100);
  assertEquals(""+shape,assertNumFound,results.numFound);
  if (assertIds != null) {
    Set<Integer> resultIds = new HashSet<>();
    for (SearchResult result : results.results) {
      resultIds.add(Integer.valueOf(result.document.get("id")));
    }
    for (int assertId : assertIds) {
      assertTrue("has " + assertId, resultIds.contains(assertId));
    }
  }
}
 

/**
 * @param doc  The doc to score
 * @return The haversine distance formula
 */
protected double distance(int doc, FunctionValues p1DV, FunctionValues p2DV) throws IOException {

  double[] p1D = new double[2];
  double[] p2D = new double[2];
  p1DV.doubleVal(doc, p1D);
  p2DV.doubleVal(doc, p2D);
  double y1;
  double x1;
  double y2;
  double x2;
  if (convertToRadians) {
    y1 = p1D[0] * DistanceUtils.DEGREES_TO_RADIANS;
    x1 = p1D[1] * DistanceUtils.DEGREES_TO_RADIANS;
    y2 = p2D[0] * DistanceUtils.DEGREES_TO_RADIANS;
    x2 = p2D[1] * DistanceUtils.DEGREES_TO_RADIANS;
  } else {
    y1 = p1D[0];
    x1 = p1D[1];
    y2 = p2D[0];
    x2 = p2D[1];
  }
  return DistanceUtils.distHaversineRAD(y1,x1,y2,x2)*radius;
}
 

protected Query makeQueryFromIndex(Shape shape) {
      // note: latitude then longitude order for LLP's methods
      if (shape instanceof Circle) {
        Circle circle = (Circle) shape;
        double radiusMeters = circle.getRadius() * DistanceUtils.DEG_TO_KM * 1000;
        return LatLonPoint.newDistanceQuery(getFieldName(),
            circle.getCenter().getY(), circle.getCenter().getX(),
            radiusMeters);
      } else if (shape instanceof Rectangle) {
        Rectangle rect = (Rectangle) shape;
        return LatLonPoint.newBoxQuery(getFieldName(),
            rect.getMinY(), rect.getMaxY(), rect.getMinX(), rect.getMaxX());
      } else if (shape instanceof Point) {
        Point point = (Point) shape;
        return LatLonPoint.newDistanceQuery(getFieldName(),
            point.getY(), point.getX(), 0);
      } else {
        throw new UnsupportedOperationException("Shape " + shape.getClass() + " is not supported by " + getClass());
      }
//      } else if (shape instanceof LucenePolygonShape) {
//        // TODO support multi-polygon
//        Polygon poly = ((LucenePolygonShape)shape).lucenePolygon;
//        return LatLonPoint.newPolygonQuery(getFieldName(), poly);
    }
 

protected Query makeQueryFromDocValues(Shape shape) {
      // note: latitude then longitude order for LLP's methods
      if (shape instanceof Circle) {
        Circle circle = (Circle) shape;
        double radiusMeters = circle.getRadius() * DistanceUtils.DEG_TO_KM * 1000;
        return LatLonDocValuesField.newSlowDistanceQuery(getFieldName(),
            circle.getCenter().getY(), circle.getCenter().getX(),
            radiusMeters);
      } else if (shape instanceof Rectangle) {
        Rectangle rect = (Rectangle) shape;
        return LatLonDocValuesField.newSlowBoxQuery(getFieldName(),
            rect.getMinY(), rect.getMaxY(), rect.getMinX(), rect.getMaxX());
      } else if (shape instanceof Point) {
        Point point = (Point) shape;
        return LatLonDocValuesField.newSlowDistanceQuery(getFieldName(),
            point.getY(), point.getX(), 0);
      } else {
        throw new UnsupportedOperationException("Shape " + shape.getClass() + " is not supported by " + getClass());
      }
//      } else if (shape instanceof LucenePolygonShape) {
//        // TODO support multi-polygon
//        Polygon poly = ((LucenePolygonShape)shape).lucenePolygon;
//        return LatLonPoint.newPolygonQuery(getFieldName(), poly);
    }
 

@Override
public Point pointOnBearing(Point from, double distDEG, double bearingDEG, SpatialContext ctx, Point reuse) {
  Geo3dPointShape geoFrom = (Geo3dPointShape) from;
  GeoPoint point = (GeoPoint) geoFrom.shape;
  double dist = DistanceUtils.DEGREES_TO_RADIANS * distDEG;
  double bearing = DistanceUtils.DEGREES_TO_RADIANS * bearingDEG;
  GeoPoint newPoint = planetModel.surfacePointOnBearing(point, dist, bearing);
  double newLat = newPoint.getLatitude() * DistanceUtils.RADIANS_TO_DEGREES;
  double newLon = newPoint.getLongitude() * DistanceUtils.RADIANS_TO_DEGREES;
  if (reuse != null) {
    reuse.reset(newLon, newLat);
    return reuse;
  }
  else {
    return ctx.getShapeFactory().pointXY(newLon, newLat);
  }
}
 

protected void initMaxLevels() {
  String mlStr = args.get(MAX_LEVELS);
  if (mlStr != null) {
    maxLevels = Integer.valueOf(mlStr);
    return;
  }

  double degrees;
  String maxDetailDistStr = args.get(MAX_DIST_ERR);
  if (maxDetailDistStr == null) {
    if (!ctx.isGeo()) {
      return;//let default to max
    }
    degrees = DistanceUtils.dist2Degrees(DEFAULT_GEO_MAX_DETAIL_KM, DistanceUtils.EARTH_MEAN_RADIUS_KM);
  } else {
    degrees = Double.parseDouble(maxDetailDistStr);
  }
  maxLevels = getLevelForDistance(degrees);
}
 

@SuppressWarnings({"unchecked"})
public SpatialScorer(LeafReaderContext readerContext, SpatialWeight w, float qWeight) throws IOException {
  super(w);
  this.weight = w;
  this.qWeight = qWeight;
  this.reader = readerContext.reader();
  this.maxDoc = reader.maxDoc();
  latVals = latSource.getValues(weight.latContext, readerContext);
  lonVals = lonSource.getValues(weight.lonContext, readerContext);

  this.lonMin = SpatialDistanceQuery.this.lonMin;
  this.lonMax = SpatialDistanceQuery.this.lonMax;
  this.lon2Min = SpatialDistanceQuery.this.lon2Min;
  this.lon2Max = SpatialDistanceQuery.this.lon2Max;
  this.latMin = SpatialDistanceQuery.this.latMin;
  this.latMax = SpatialDistanceQuery.this.latMax;
  this.lon2 = SpatialDistanceQuery.this.lon2;
  this.calcDist = SpatialDistanceQuery.this.calcDist;

  this.latCenterRad = SpatialDistanceQuery.this.latCenter * DistanceUtils.DEGREES_TO_RADIANS;
  this.lonCenterRad = SpatialDistanceQuery.this.lonCenter * DistanceUtils.DEGREES_TO_RADIANS;
  this.latCenterRad_cos = this.calcDist ? Math.cos(latCenterRad) : 0;
  this.dist = SpatialDistanceQuery.this.dist;
  this.planetRadius = SpatialDistanceQuery.this.planetRadius;

}
 

double dist(double lat, double lon) {
  double latRad = lat * DistanceUtils.DEGREES_TO_RADIANS;
  double lonRad = lon * DistanceUtils.DEGREES_TO_RADIANS;
  
  // haversine, specialized to avoid a cos() call on latCenterRad
  double diffX = latCenterRad - latRad;
  double diffY = lonCenterRad - lonRad;
  double hsinX = Math.sin(diffX * 0.5);
  double hsinY = Math.sin(diffY * 0.5);
  double h = hsinX * hsinX +
          (latCenterRad_cos * Math.cos(latRad) * hsinY * hsinY);
  double result = (planetRadius * 2 * Math.atan2(Math.sqrt(h), Math.sqrt(1 - h)));

  // save the results of this calculation
  lastDistDoc = doc;
  lastDist = result;
  
  return result;
}
 

/**
 * Convert a value to degrees
 *
 * @param v     value
 * @param units UoM IRI of the unit
 * @return degrees as a double
 * @throws ValueExprEvaluationException
 */
public static double convertToDegrees(double v, IRI units) throws ValueExprEvaluationException {
	double degs;

	if (GEOF.UOM_DEGREE.equals(units)) {
		degs = v;
	} else if (GEOF.UOM_RADIAN.equals(units)) {
		degs = DistanceUtils.toDegrees(v);
	} else if (GEOF.UOM_UNITY.equals(units)) {
		degs = v * 180.0;
	} else if (GEOF.UOM_METRE.equals(units)) {
		degs = DistanceUtils.dist2Degrees(v / 1000.0, DistanceUtils.EARTH_MEAN_RADIUS_KM);
	} else {
		throw new ValueExprEvaluationException("Invalid unit of measurement: " + units);
	}
	return degs;
}
 

@Override
protected Shape randomIndexedShape() {
  Rectangle world = ctx.getWorldBounds();
  if (random().nextInt(10) == 0) // increased chance of getting one of these
    return world;

  int worldWidth = (int) Math.round(world.getWidth());
  int deltaLeft = nextIntInclusive(worldWidth);
  int deltaRight = nextIntInclusive(worldWidth - deltaLeft);
  int worldHeight = (int) Math.round(world.getHeight());
  int deltaTop = nextIntInclusive(worldHeight);
  int deltaBottom = nextIntInclusive(worldHeight - deltaTop);
  if (ctx.isGeo() && (deltaLeft != 0 || deltaRight != 0)) {
    //if geo & doesn't world-wrap, we shift randomly to potentially cross dateline
    int shift = nextIntInclusive(360);
    return ctx.getShapeFactory().rect(
        DistanceUtils.normLonDEG(world.getMinX() + deltaLeft + shift),
        DistanceUtils.normLonDEG(world.getMaxX() - deltaRight + shift),
        world.getMinY() + deltaBottom, world.getMaxY() - deltaTop);
  } else {
    return ctx.getShapeFactory().rect(
        world.getMinX() + deltaLeft, world.getMaxX() - deltaRight,
        world.getMinY() + deltaBottom, world.getMaxY() - deltaTop);
  }
}
 

@Override
public double getDistance() {
	String geohash = (String) ((ElasticsearchDocument) getDocument()).getSource().get(geoPointField);
	GeoPoint dstPoint = GeoPoint.fromGeohash(geohash);

	double unitDist = GeoDistance.ARC.calculate(srcPoint.getLat(), srcPoint.getLon(), dstPoint.getLat(),
			dstPoint.getLon(), unit);
	double distance;
	if (GEOF.UOM_METRE.equals(units)) {
		distance = unit.toMeters(unitDist);
	} else if (GEOF.UOM_DEGREE.equals(units)) {
		distance = unitDist / unit.getDistancePerDegree();
	} else if (GEOF.UOM_RADIAN.equals(units)) {
		distance = DistanceUtils.dist2Radians(unit.convert(unitDist, DistanceUnit.KILOMETERS),
				DistanceUtils.EARTH_MEAN_RADIUS_KM);
	} else if (GEOF.UOM_UNITY.equals(units)) {
		distance = unit.convert(unitDist, DistanceUnit.KILOMETERS) / (Math.PI * DistanceUtils.EARTH_MEAN_RADIUS_KM);
	} else {
		throw new UnsupportedOperationException("Unsupported units: " + units);
	}
	return distance;
}
 

private void validateHeatmapResult(Rectangle inputRange, int facetLevel, HeatmapFacetCounter.Heatmap heatmap)
    throws IOException {
  final Rectangle heatRect = heatmap.region;
  assertTrue(heatRect.relate(inputRange) == SpatialRelation.CONTAINS || heatRect.equals(inputRange));
  final double cellWidth = heatRect.getWidth() / heatmap.columns;
  final double cellHeight = heatRect.getHeight() / heatmap.rows;
  for (int c = 0; c < heatmap.columns; c++) {
    for (int r = 0; r < heatmap.rows; r++) {
      final int facetCount = heatmap.getCount(c, r);
      double x = DistanceUtils.normLonDEG(heatRect.getMinX() + c * cellWidth + cellWidth / 2);
      double y = DistanceUtils.normLatDEG(heatRect.getMinY() + r * cellHeight + cellHeight / 2);
      Point pt =  shapeFactory.pointXY(x, y);
      assertEquals(countMatchingDocsAtLevel(pt, facetLevel), facetCount);
    }
  }
}
 

private ShapeBuilder getCircleBuilder(GeoCircle geoCircle) {
    // NOTE: as of ES7, storing circles is no longer supported so we need approximate the circle with a polygon
    double radius = geoCircle.getRadius();
    double maxSideLengthKm = getSearchIndex().getConfig().getGeocircleToPolygonSideLength();
    maxSideLengthKm = Math.min(radius, maxSideLengthKm);

    // calculate how many points we need to use given the length of a polygon side
    int numberOfPoints = (int) Math.ceil(Math.PI / Math.asin((maxSideLengthKm / (2 * radius))));
    numberOfPoints = Math.min(numberOfPoints, getSearchIndex().getConfig().getGeocircleToPolygonMaxNumSides());

    // Given the number of sides, loop through slices of 360 degrees and calculate the lat/lon at that radius and heading
    SpatialContext spatialContext = SpatialContext.GEO;
    DistanceCalculator distanceCalculator = spatialContext.getDistCalc();
    Point centerPoint = spatialContext.getShapeFactory().pointXY(DistanceUtils.normLonDEG(geoCircle.getLongitude()), DistanceUtils.normLatDEG(geoCircle.getLatitude()));
    ArrayList<GeoPoint> points = new ArrayList<>();
    for (float angle = 360; angle > 0; angle -= 360.0 / numberOfPoints) {
        Point point = distanceCalculator.pointOnBearing(centerPoint, geoCircle.getRadius() * KM_TO_DEG, angle, spatialContext, null);
        points.add(new GeoPoint(point.getY(), point.getX()));
    }

    // Polygons must start/end at the same point, so add the first point onto the end
    points.add(points.get(0));

    return getPolygonBuilder(new GeoPolygon(points, geoCircle.getDescription()));
}
 

@Override
public void reset(double x, double y, double radiusDEG) {
  shape = GeoCircleFactory.makeGeoCircle(shape.getPlanetModel(),
      y * DistanceUtils.DEGREES_TO_RADIANS,
      x * DistanceUtils.DEGREES_TO_RADIANS,
      radiusDEG * DistanceUtils.DEGREES_TO_RADIANS);
  center = null;
  boundingBox = null;
}
 

@Test
public void testPrecision() throws IOException{
  init(GeohashPrefixTree.getMaxLevelsPossible());

  Point iPt = ctx.getShapeFactory().pointXY(2.8028712999999925, 48.3708044);//lon, lat
  addDocument(newDoc("iPt", iPt));
  commit();

  Point qPt = ctx.getShapeFactory().pointXY(2.4632387000000335, 48.6003516);

  final double KM2DEG = DistanceUtils.dist2Degrees(1, DistanceUtils.EARTH_MEAN_RADIUS_KM);
  final double DEG2KM = 1 / KM2DEG;

  final double DIST = 35.75;//35.7499...
  assertEquals(DIST, ctx.getDistCalc().distance(iPt, qPt) * DEG2KM, 0.001);

  //distErrPct will affect the query shape precision. The indexed precision
  // was set to nearly zilch via init(GeohashPrefixTree.getMaxLevelsPossible());
  final double distErrPct = 0.025; //the suggested default, by the way
  final double distMult = 1+distErrPct;

  assertTrue(35.74*distMult >= DIST);
  checkHits(q(qPt, 35.74 * KM2DEG, distErrPct), 1, null);

  assertTrue(30*distMult < DIST);
  checkHits(q(qPt, 30 * KM2DEG, distErrPct), 0, null);

  assertTrue(33*distMult < DIST);
  checkHits(q(qPt, 33 * KM2DEG, distErrPct), 0, null);

  assertTrue(34*distMult < DIST);
  checkHits(q(qPt, 34 * KM2DEG, distErrPct), 0, null);
}
 

@Test
public void testOneMeterPrecision() {
  init(GeohashPrefixTree.getMaxLevelsPossible());
  GeohashPrefixTree grid = (GeohashPrefixTree) ((RecursivePrefixTreeStrategy) strategy).getGrid();
  //DWS: I know this to be true.  11 is needed for one meter
  double degrees = DistanceUtils.dist2Degrees(0.001, DistanceUtils.EARTH_MEAN_RADIUS_KM);
  assertEquals(11, grid.getLevelForDistance(degrees));
}
 

protected Rectangle randomRectangle(Rectangle bounds) {
  double[] xNewStartAndWidth = randomSubRange(bounds.getMinX(), bounds.getWidth());
  double xMin = xNewStartAndWidth[0];
  double xMax = xMin + xNewStartAndWidth[1];
  if (bounds.getCrossesDateLine()) {
    xMin = DistanceUtils.normLonDEG(xMin);
    xMax = DistanceUtils.normLonDEG(xMax);
  }

  double[] yNewStartAndHeight = randomSubRange(bounds.getMinY(), bounds.getHeight());
  double yMin = yNewStartAndHeight[0];
  double yMax = yMin + yNewStartAndHeight[1];

  return ctx.getShapeFactory().rect(xMin, xMax, yMin, yMax);
}
 

public static final double toMiles(double distance, IRI units) {
	final double miles;
	if (GEOF.UOM_METRE.equals(units)) {
		miles = DistanceUtils.KM_TO_MILES * distance / 1000.0;
	} else if (GEOF.UOM_DEGREE.equals(units)) {
		miles = DistanceUtils.degrees2Dist(distance, DistanceUtils.EARTH_MEAN_RADIUS_MI);
	} else if (GEOF.UOM_RADIAN.equals(units)) {
		miles = DistanceUtils.radians2Dist(distance, DistanceUtils.EARTH_MEAN_RADIUS_MI);
	} else if (GEOF.UOM_UNITY.equals(units)) {
		miles = distance * Math.PI * DistanceUtils.EARTH_MEAN_RADIUS_MI;
	} else {
		throw new IllegalArgumentException("Unsupported units: " + units);
	}
	return miles;
}
 

@Override
public DoubleValuesSource makeDistanceValueSource(Point queryPoint, double multiplier) {
  if (!docValues) {
    throw new SolrException(SolrException.ErrorCode.BAD_REQUEST,
        getFieldName() + " must have docValues enabled to support this feature");
  }
  // Internally, the distance from LatLonPointSortField/Comparator is in meters. So we must also go from meters to
  //  degrees, which is what Lucene spatial-extras is oriented around.
  return new DistanceSortValueSource(getFieldName(), queryPoint,
      DistanceUtils.KM_TO_DEG / 1000.0 * multiplier);
}
 

/**
 * Implemented for compatibility with geofilt & bbox query parsers:
 * {@link SpatialQueryable}.
 */
@Override
public Query createSpatialQuery(QParser parser, SpatialOptions options) {
  Point pt = SpatialUtils.parsePointSolrException(options.pointStr, ctx);

  double distDeg = DistanceUtils.dist2Degrees(options.distance, options.radius);

  Shape shape = ctx.makeCircle(pt, distDeg);
  if (options.bbox)
    shape = shape.getBoundingBox();

  SpatialArgs spatialArgs = new SpatialArgs(SpatialOperation.Intersects, shape);
  return getQueryFromSpatialArgs(parser, options.field, spatialArgs);
}
 

@Test
public void testHaversine() throws Exception {
  clearIndex();
  assertU(adoc("id", "1", "x_td", "0", "y_td", "0", "gh_s1", GeohashUtils.encodeLatLon(32.7693246, -79.9289094)));
  assertU(adoc("id", "2", "x_td", "0", "y_td", String.valueOf(Math.PI / 2), "gh_s1", GeohashUtils.encodeLatLon(32.7693246, -78.9289094)));
  assertU(adoc("id", "3", "x_td", String.valueOf(Math.PI / 2), "y_td", String.valueOf(Math.PI / 2), "gh_s1", GeohashUtils.encodeLatLon(32.7693246, -80.9289094)));
  assertU(adoc("id", "4", "x_td", String.valueOf(Math.PI / 4), "y_td", String.valueOf(Math.PI / 4), "gh_s1", GeohashUtils.encodeLatLon(32.7693246, -81.9289094)));
  assertU(adoc("id", "5", "x_td", "45.0", "y_td", "45.0",
          "gh_s1", GeohashUtils.encodeLatLon(32.7693246, -81.9289094)));
  assertU(adoc("id", "6", "point_hash", "32.5, -79.0", "point", "32.5, -79.0"));
  assertU(adoc("id", "7", "point_hash", "32.6, -78.0", "point", "32.6, -78.0"));
  assertU(commit());
  //Get the haversine distance between the point 0,0 and the docs above assuming a radius of 1
  assertQ(req("fl", "*,score", "q", "{!func}hsin(1, false, x_td, y_td, 0, 0)", "fq", "id:1"), "//float[@name='score']='0.0'");
  assertQ(req("fl", "*,score", "q", "{!func}hsin(1, false, x_td, y_td, 0, 0)", "fq", "id:2"), "//float[@name='score']='" + (float) (Math.PI / 2) + "'");
  assertQ(req("fl", "*,score", "q", "{!func}hsin(1, false, x_td, y_td, 0, 0)", "fq", "id:3"), "//float[@name='score']='" + (float) (Math.PI / 2) + "'");
  assertQ(req("fl", "*,score", "q", "{!func}hsin(1, false, x_td, y_td, 0, 0)", "fq", "id:4"), "//float[@name='score']='1.0471976'");
  assertQ(req("fl", "*,score", "q", "{!func}hsin(1, true, x_td, y_td, 0, 0)", "fq", "id:5"), "//float[@name='score']='1.0471976'");
  //SOLR-2114
  assertQ(req("fl", "*,score", "q", "{!func}hsin(6371.009, true, point, vector(0, 0))", "fq", "id:6"), "//float[@name='score']='8977.814'");
  
  //Geo Hash Haversine
  //Can verify here: http://www.movable-type.co.uk/scripts/latlong.html, but they use a slightly different radius for the earth, so just be close
  //note: using assertJQ because it supports numeric deltas, and by default too
  assertJQ(req("fl", "*,score", "q", "{!func}ghhsin(" + DistanceUtils.EARTH_MEAN_RADIUS_KM + ", gh_s1, \"" + GeohashUtils.encodeLatLon(32, -79) + "\",)", "fq", "id:1"),
      "/response/docs/[0]/score==122.171875");

  assertQ(req("fl", "id,point_hash,score", "q", "{!func}recip(ghhsin(" + DistanceUtils.EARTH_MEAN_RADIUS_KM + ", point_hash, \"" + GeohashUtils.encodeLatLon(32, -79) + "\"), 1, 1, 0)"),
          "//*[@numFound='7']", 
          "//result/doc[1]/str[@name='id'][.='6']",
          "//result/doc[2]/str[@name='id'][.='7']"//all the rest don't matter
          );


  assertJQ(req("fl", "*,score", "q", "{!func}ghhsin(" + DistanceUtils.EARTH_MEAN_RADIUS_KM + ", gh_s1, geohash(32, -79))", "fq", "id:1"),
      "/response/docs/[0]/score==122.171875");

}
 

@Test @Repeat(iterations = 10)
public void testLLPDecodeIsStableAndPrecise() throws Exception {
  // test that LatLonPointSpatialField decode of docValue will round-trip (re-index then re-decode) to the same value
  @SuppressWarnings({"resource", "IOResourceOpenedButNotSafelyClosed"})
  SolrClient client = new EmbeddedSolrServer(h.getCore());// do NOT close it; it will close Solr

  final String fld = "llp_1_dv_dvasst";
  String ptOrig = GeoTestUtil.nextLatitude() + "," + GeoTestUtil.nextLongitude();
  assertU(adoc("id", "0", fld, ptOrig));
  assertU(commit());
  // retrieve it (probably less precision)
  String ptDecoded1 = (String) client.query(params("q", "id:0")).getResults().get(0).get(fld);
  // now write it back
  assertU(adoc("id", "0", fld, ptDecoded1));
  assertU(commit());
  // retrieve it; assert that it's the same as written
  String ptDecoded2 = (String) client.query(params("q", "id:0")).getResults().get(0).get(fld);
  assertEquals("orig:" + ptOrig, ptDecoded1, ptDecoded2);

  // test that the representation is pretty accurate
  final Point ptOrigObj = SpatialUtils.parsePoint(ptOrig, SpatialContext.GEO);
  final Point ptDecodedObj = SpatialUtils.parsePoint(ptDecoded1, SpatialContext.GEO);
  double deltaCentimeters = SpatialContext.GEO.calcDistance(ptOrigObj, ptDecodedObj) * DistanceUtils.DEG_TO_KM * 1000.0 * 100.0;
  //See javadocs of LatLonDocValuesField for these constants
  final Point absErrorPt = SpatialContext.GEO.getShapeFactory().pointXY(8.381903171539307E-8, 4.190951585769653E-8);
  double deltaCentimetersMax
      = SpatialContext.GEO.calcDistance(absErrorPt, 0,0) * DistanceUtils.DEG_TO_KM * 1000.0 * 100.0;
  assertEquals(1.0420371840922256, deltaCentimetersMax, 0.0);// just so that we see it in black & white in the test

  //max found by trial & error.  If we used 8 decimal places then we could get down to 1.04cm accuracy but then we
  // lose the ability to round-trip -- 40 would become 39.99999997  (ugh).
  assertTrue("deltaCm too high: " + deltaCentimeters, deltaCentimeters < 1.41);
  // Pt(x=105.29894270124083,y=-0.4371673760042398) to  Pt(x=105.2989428,y=-0.4371673) is 1.38568
}
 

@Test
public void testIntersectFilter() throws Exception {
  setupDocs();
  
  //Try some edge cases
  checkHits(fieldName, "1,1", 175, DistanceUtils.EARTH_MEAN_RADIUS_KM, 3, 5, 6, 7);
  checkHits(fieldName, "0,179.8", 200, DistanceUtils.EARTH_MEAN_RADIUS_KM, 2, 8, 9);
  checkHits(fieldName, "89.8, 50", 200, DistanceUtils.EARTH_MEAN_RADIUS_KM, 2, 10, 11);//this goes over the north pole
  checkHits(fieldName, "-89.8, 50", 200, DistanceUtils.EARTH_MEAN_RADIUS_KM, 2, 12, 13);//this goes over the south pole
  //try some normal cases
  checkHits(fieldName, "33.0,-80.0", 300, DistanceUtils.EARTH_MEAN_RADIUS_KM, 2);
  //large distance
  checkHits(fieldName, "1,1", 5000, DistanceUtils.EARTH_MEAN_RADIUS_KM, 3, 5, 6, 7);
  //Because we are generating a box based on the west/east longitudes and the south/north latitudes, which then
  //translates to a range query, which is slightly more inclusive.  Thus, even though 0.0 is 15.725 kms away,
  //it will be included, b/c of the box calculation.
  checkHits(fieldName, false, "0.1,0.1", 15, DistanceUtils.EARTH_MEAN_RADIUS_KM, 2, 5, 6);

  //try some more
  clearIndex();
  assertU(adoc("id", "14", fieldName, "0,5"));
  assertU(adoc("id", "15", fieldName, "0,15"));
  //3000KM from 0,0, see http://www.movable-type.co.uk/scripts/latlong.html
  assertU(adoc("id", "16", fieldName, "18.71111,19.79750"));
  assertU(adoc("id", "17", fieldName, "44.043900,-95.436643"));
  assertU(commit());

  checkHits(fieldName, "0,0", 1000, DistanceUtils.EARTH_MEAN_RADIUS_KM, 1, 14);
  checkHits(fieldName, "0,0", 2000, DistanceUtils.EARTH_MEAN_RADIUS_KM, 2, 14, 15);
  checkHits(fieldName, false, "0,0", 3000, DistanceUtils.EARTH_MEAN_RADIUS_KM, 3, 14, 15, 16);
  checkHits(fieldName, "0,0", 3001, DistanceUtils.EARTH_MEAN_RADIUS_KM, 3, 14, 15, 16);
  checkHits(fieldName, "0,0", 3000.1, DistanceUtils.EARTH_MEAN_RADIUS_KM, 3, 14, 15, 16);

  //really fine grained distance and reflects some of the vagaries of how we are calculating the box
  checkHits(fieldName, "43.517030,-96.789603", 109, DistanceUtils.EARTH_MEAN_RADIUS_KM, 0);

  //falls outside of the real distance, but inside the bounding box
  checkHits(fieldName, true,  "43.517030,-96.789603", 110, DistanceUtils.EARTH_MEAN_RADIUS_KM, 0);
  checkHits(fieldName, false, "43.517030,-96.789603", 110, DistanceUtils.EARTH_MEAN_RADIUS_KM, 1, 17);
}
 

/**
 * Return a lucene query that finds all points within the given radius from the given point
 */
public static Query findWithin(double longitude, double latitude, double radiusMiles) {
  // Covert the radius in miles to a radius in degrees
  double radiusDEG = DistanceUtils.dist2Degrees(radiusMiles, EARTH_MEAN_RADIUS_MI);

  // Create a query that looks for all points within a circle around the given point
  SpatialArgs args = new SpatialArgs(SpatialOperation.IsWithin,
      new GeoCircle(createPoint(longitude, latitude), radiusDEG, CONTEXT));
  return STRATEGY.makeQuery(args);
}
 

private void computeDistance(Bucket bucket, Bucket potentialNeighbor, List<Bucket> revisit, ReduceContext reduceContext) {
    if (potentialNeighbor.visited) {
        return;
    }

    double neighborDistance = GeoUtils.arcDistance(
            bucket.centroid.lat(),
            bucket.centroid.lon(),
            potentialNeighbor.centroid.lat(),
            potentialNeighbor.centroid.lon()
    );

    double avgLat = (bucket.centroid.lat() + potentialNeighbor.centroid.lat()) / 2;

    double fixedRadius = radius * Math.cos(DistanceUtils.toRadians(avgLat));

    if (neighborDistance <= fixedRadius) {
        potentialNeighbor.visited = true;
        long mergedDocCount = bucket.docCount + potentialNeighbor.docCount;
        double newCentroidLat = (bucket.centroid.getLat() * bucket.docCount +
                potentialNeighbor.centroid.getLat() * potentialNeighbor.docCount) / mergedDocCount;
        double newCentroidLon = (bucket.centroid.getLon() * bucket.docCount +
                potentialNeighbor.centroid.getLon() * potentialNeighbor.docCount) / mergedDocCount;
        bucket.centroid = new GeoPoint(newCentroidLat, newCentroidLon);

        bucket.docCount = mergedDocCount;
        List<InternalAggregations> aggregationsList = new ArrayList<>();
        aggregationsList.add(bucket.aggregations);
        aggregationsList.add(potentialNeighbor.aggregations);
        bucket.aggregations = InternalAggregations.reduce(aggregationsList, reduceContext);
        bucket.geohashesList.add(potentialNeighbor.geohashAsLong);
    } else if (revisit != null && ratio > 0 && neighborDistance / fixedRadius < ratio ) {
        revisit.add(potentialNeighbor);
    }
}
 

@Override
public double distance(Point from, double toX, double toY) {
  GeoPoint fromGeoPoint;
  if (from instanceof Geo3dPointShape) {
    fromGeoPoint = (((Geo3dPointShape) from).shape).getCenter();
  } else {
    fromGeoPoint = new GeoPoint(planetModel,
        from.getY() * DistanceUtils.DEGREES_TO_RADIANS,
        from.getX() * DistanceUtils.DEGREES_TO_RADIANS);
  }
  GeoPoint toGeoPoint = new GeoPoint(planetModel,
      toY * DistanceUtils.DEGREES_TO_RADIANS,
      toX * DistanceUtils.DEGREES_TO_RADIANS);
  return planetModel.surfaceDistance(fromGeoPoint, toGeoPoint) * DistanceUtils.RADIANS_TO_DEGREES;
}
 

@Override
public double getDistanceForLevel(int level) {
    if (level == 0) {
        return 180;
    }
    return S2Projections.MAX_WIDTH.getValue(arity * (level - 1)) * DistanceUtils.RADIANS_TO_DEGREES;
}
 

@Override
public int getLevelForDistance(double dist) {
    if (dist == 0){
        return maxLevels;
    }
    int level =  S2Projections.MAX_WIDTH.getMinLevel(dist * DistanceUtils.DEGREES_TO_RADIANS);
    int roundLevel = level % arity != 0 ? 1 : 0;
    level = level/arity + roundLevel;
    return Math.min(maxLevels, level + 1);
}
 

@Override
protected Iterable<? extends DocumentDistance> geoQuery(final IRI geoProperty, Point p, final IRI units,
		double distance, String distanceVar, Var contextVar) throws MalformedQueryException, IOException {
	double unitDist;
	final DistanceUnit unit;
	if (GEOF.UOM_METRE.equals(units)) {
		unit = DistanceUnit.METERS;
		unitDist = distance;
	} else if (GEOF.UOM_DEGREE.equals(units)) {
		unit = DistanceUnit.KILOMETERS;
		unitDist = unit.getDistancePerDegree() * distance;
	} else if (GEOF.UOM_RADIAN.equals(units)) {
		unit = DistanceUnit.KILOMETERS;
		unitDist = DistanceUtils.radians2Dist(distance, DistanceUtils.EARTH_MEAN_RADIUS_KM);
	} else if (GEOF.UOM_UNITY.equals(units)) {
		unit = DistanceUnit.KILOMETERS;
		unitDist = distance * Math.PI * DistanceUtils.EARTH_MEAN_RADIUS_KM;
	} else {
		throw new MalformedQueryException("Unsupported units: " + units);
	}

	double lat = p.getY();
	double lon = p.getX();
	final String fieldName = toGeoPointFieldName(SearchFields.getPropertyField(geoProperty));
	QueryBuilder qb = QueryBuilders.functionScoreQuery(
			QueryBuilders.geoDistanceQuery(fieldName).point(lat, lon).distance(unitDist, unit),
			ScoreFunctionBuilders.linearDecayFunction(fieldName, GeohashUtils.encodeLatLon(lat, lon),
					new DistanceUnit.Distance(unitDist, unit).toString()));
	if (contextVar != null) {
		qb = addContextTerm(qb, (Resource) contextVar.getValue());
	}

	SearchRequestBuilder request = client.prepareSearch();
	SearchHits hits = search(request, qb);
	final GeoPoint srcPoint = new GeoPoint(lat, lon);
	return Iterables.transform(hits, (Function<SearchHit, DocumentDistance>) hit -> {
		return new ElasticsearchDocumentDistance(hit, geoContextMapper, fieldName, units, srcPoint, unit);
	});
}