[GIS] How to calculate whether a point is in a shape using geodetic calculations

geotoolspoint-in-polygon

I’m new to GIS so forgive me if my terminology isn’t completely correct. Really liking the options available in the GeoTools library but I’ve hit a stumbling block that I’m hoping someone will be able to advise on. I’ve looked at the classes in the library but can’t find a method that is giving me what I need. I have a scenario where I have a point and a polygon defined by latitude and longitude (WGS84). I need to determine whether the point is in (or on the boundary) of the polygon.

I have tried using the Intersects method but it looks to be treating the polygon as if it is on a 2D plane. The bottom left and right corners of my polygon have the same latitude value but they are spaced longitudinally by 10 degrees, so it’s quite a long edge. There are no intermediate points along this edge of the polygon. The Intersects method appears to use a straight line between these two points, effectively taking this bottom edge under ground when on a globe.

What I need is for this polygon’s edges to be interpreted as great circle lines between points on the surface of the globe, thereby creating slight curves, so when we determine whether a point is in the polygon, we are using the shape with slightly curved edges. Obviously the greater distance between polygon points the greater the curve.

Is there a method to achieve this in the GeoTools library? If not, does anyone have alternative open source Java libraries that can solve what must be a common request?

Best Answer

You need the JTS Densifier which takes a geometry and returns the same line but with more points in it.

    GeometryFactory gf = new GeometryFactory();
    LineString l = gf.createLineString(new Coordinate[] {new Coordinate(-10, 50),new Coordinate(10, 50)});
    WKTWriter2 writer = new WKTWriter2();
    System.out.println(writer.writeFormatted(l));
    LineString dl = (LineString) Densifier.densify(l, 1);
    System.out.println(writer.writeFormatted(dl));

Which gives you the following:

LINESTRING (-10 50, 10 50)
LINESTRING (-10 50, -9.047619047619047 50, -8.095238095238095 50, -7.142857142857143 50, -6.190476190476191 50, -5.238095238095238 50, -4.2857142857142865 50, -3.333333333333334 50, -2.3809523809523814 50, -1.4285714285714288 50, -0.4761904761904763 50, 0.4761904761904763 50, 1.428571428571427 50, 2.3809523809523796 50, 3.333333333333332 50, 4.285714285714285 50, 5.238095238095237 50, 6.19047619047619 50, 7.142857142857142 50, 8.095238095238095 50, 9.047619047619047 50, 10 50)

How these appear on your map will depend on your projection, if you are using a plate carree projection (e.g 4326) then they are identical if however you use a Transverse Mercator like OSGB (27700) the 2nd one is curved.

Related Solutions

[GIS] I nead to find the nearest places based on latitude/longitude points

Yes, dealing with distance (or area) in unprojected latitude/longitude coordinates is fraught with peril. The length of a degree of latitude is always the same (60 nautical miles), but the length of a degree of longitude varies with the cosine of the latitude.

To do it properly, you want to project your coordinates into some plane coordinate system and do your comparisons in that. (Which plane coordinate system should you use? Not a simple answer, depends where on earth your stuff is.) Geotools must have a library for that, or since you're working in java, check out proj4j.

Alternatively, it's kind of a hack, but if all you are doing is calculating straight-line distances, you could create sort of pseudo-coordinates in meters where y=latitude*60*cf (where cf is the conversion factor to get from nautical miles to meters) and x=longitude*60*cf/cos(latitude). You're assuming a spherical earth, which it isn't. And it will severely break down if your distances are very big. But for 500m it ought to work. If you do it this way it might just be easier to calculate distance yourself, using d=sqr((x2-x1)^2 + (y2-y1)^2)).

Better would be to use Vincenty's formula, which doesn't assume a sphere. Poke around for that, you may find a library that will do it.

[GIS] Finding distance in meters between point and path in lat/lon using GeoTools

The code from the previous question does in fact solve this problem once you have explored the issue a bit.

Basically your line doesn't have enough points so it "cuts" the corner across the more localised projection that is used to allow JTS to proceed as if the Earth is flat. This needs one additional line of code:

route = (LineString) Densifier.densify(route, .01);

The brown line is your original 3 point line and the blue line is a densified line (that has points added so it projects properly), shown in a local projection (EPSG:2470). As you can see the distance to the brown line to the point is much shorter than it should be.

So here is some suitably modified code:

public double getDistance(Point point, LineString line) {
double dist = -1.0;
try {
  String code = "AUTO:42001," + point.getX() + "," + point.getY();
  CoordinateReferenceSystem auto = CRS.decode(code);
  // auto = CRS.decode("epsg:2470");
  MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84, auto);
  Geometry g3 = JTS.transform(line, transform);
  Geometry g4 = JTS.transform(point, transform);

  Coordinate[] c = DistanceOp.nearestPoints(g4, g3);

  Coordinate c1 = new Coordinate();
  //System.out.println(c[1].distance(g4.getCoordinate()));
  JTS.transform(c[1], c1, transform.inverse());
  //System.out.println(geometryFactory.createPoint(c1));
  dist = JTS.orthodromicDistance(point.getCoordinate(), c1, DefaultGeographicCRS.WGS84);
} catch (Exception e) {
  e.printStackTrace();
}
return dist;

}

And incase you are wondering about the difference between the "flat distance" (c[1].distance(g4.getCoordinate())) is 760.731km and the "curved distance" (JTS.orthodromicDistance(point.getCoordinate(), c1, DefaultGeographicCRS.WGS84);) is 760.980km. For reference the distance tool in QGIS gives 760.919 km.

Best Answer

Related Solutions

[GIS] I nead to find the nearest places based on latitude/longitude points

[GIS] Finding distance in meters between point and path in lat/lon using GeoTools

Related Question