1

New to GeoTools and GIS and I am trying to calculate distance between Mumbai and Durban using GeoTools library. I am getting close to accurate results for small distances but when i go for bigger ones,the calculation is way too offcourse by 2000 km, i dont completely understand the CRS system .Below is my Code to calculate the distance between Mumbai and Durban

    Coordinate source = new Coordinate(19.0760, 72.8777);   ///Mumbai Lat Long
    Coordinate destination1 = new Coordinate(-29.883333, 31.049999); //Durban Lat Long

    GeometryFactory geometryFactory = new GeometryFactory();
    Geometry point1 = geometryFactory.createPoint(source);
    Geometry point2 = geometryFactory.createPoint(destination1);

    CoordinateReferenceSystem auto = auto = CRS.decode("AUTO:42001,13.45,52.3");
    MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84, auto);

    Geometry g3 = JTS.transform(point1, transform);
    Geometry g4 = JTS.transform(point2, transform);

    double distance = g3.distance(g4);

1 Answers1

8

This is what happens when you copy code blindly from stackexchange questions without reading the question it was based on which explains why.

All the times I've answered that question (and posted code like that) the questioner is trying to use lat/lon coordinates in degrees to measure a short distance in metres. The trick shown in your question creates an automatic UTM projection centred on the position specified after the "AUTO:42001," bit (in your case 52N 13E) - this needs to be the centre of the area you are interested in, so in your case those values are probably wrong anyway.

But you aren't interested in a small region Mumbai to Durban is a significant way around the Earth so you need to allow for the curvature of the Earth's surface. Also you aren't trying to do something difficult for which JTS is the only source of process (e.g buffering). In this case you should use the GeodeticCalculator which takes the shape of the Earth into account using the library from C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43–55 (2013).

Anyway enough explanation that no one will read in the future, here's the code:

  public static void main(String[] args) {
    DefaultGeographicCRS crs = DefaultGeographicCRS.WGS84;
    if (args.length != 4) {
      System.err.println("Need 4 numbers lat_1 lon_1 lat_2 lon_2");
      return;
    }
    GeometryFactory geomFactory = new GeometryFactory();
    Point[] points = new Point[2];
    for (int i = 0, k = 0; i < 2; i++, k += 2) {
      double x = Double.valueOf(args[k]);
      double y = Double.valueOf(args[k + 1]);
      if (CRS.getAxisOrder(crs).equals(AxisOrder.NORTH_EAST)) {
        System.out.println("working with a lat/lon crs");
        points[i] = geomFactory.createPoint(new Coordinate(x, y));
      } else {
        System.out.println("working with a lon/lat crs");
        points[i] = geomFactory.createPoint(new Coordinate(y, x));
      }

    }

    double distance = 0.0;

    GeodeticCalculator calc = new GeodeticCalculator(crs);
    calc.setStartingGeographicPoint(points[0].getX(), points[0].getY());
    calc.setDestinationGeographicPoint(points[1].getX(), points[1].getY());

    distance = calc.getOrthodromicDistance();
    double bearing = calc.getAzimuth();

    Quantity<Length> dist = Quantities.getQuantity(distance, SI.METRE);
    System.out.println(dist.to(MetricPrefix.KILO(SI.METRE)).getValue() + " Km");
    System.out.println(dist.to(USCustomary.MILE).getValue() + " miles");
    System.out.println("Bearing " + bearing + " degrees");
  }

Giving:

working with a lon/lat crs
POINT (72.8777 19.076)
POINT (31.049999 -29.883333)
7032.866960793305 Km
4370.020928274692 miles
Bearing -139.53428618565218 degrees
Ian Turton
  • 10,018
  • 1
  • 28
  • 47
  • This is exactly what i wanted !!!. Thank you Ian. appreciated – Nikhil Karanjkar Feb 22 '20 at 15:09
  • i also wanted to find out the intersecting distance between a linestring and polygon. The way i am trying to do this is by projecting the geometries to auto projections from wGs84 by using a centroid of the polygon (but i think i should take a centroid of both geometries and i am currently not sure how to do that), and then getting the distance of the intersecting geometry. But i found that it is off course. Can you advice me how should i approach on this ? – Nikhil Karanjkar Feb 22 '20 at 17:29
  • @Ian , what libraries are you using for the code ? " Quantity dist = Quantities.getQuantity(distance, SI.METRE);" – Jeryl Cook Jul 27 '21 at 09:57
  • Is this still valid for small distances like few meters ? – Magno C Mar 27 '23 at 03:21
  • I think it's supposed to work down to millimetres – Ian Turton Mar 27 '23 at 07:32
  • Below the libraries that I used: import org.locationtech.jts.geom.GeometryFactory; import org.geotools.referencing.CRS; import org.geotools.referencing.GeodeticCalculator; import org.geotools.referencing.CRS.AxisOrder; import org.geotools.referencing.crs.DefaultGeographicCRS; import si.uom.SI; import tech.units.indriya.quantity.Quantities; import javax.measure.MetricPrefix; import javax.measure.Quantity; import javax.measure.quantity.Length; – Marcos Roberto Silva Apr 16 '23 at 11:30