Saturday 10 June 2017

Using Linestring distance in meters with GeoTools?


I want to use GeoTools for working with gis data. I'm most concerned about distance calculations, not so much about drawing a map. What I still couldn't figure out is how to make distance calculations work smoothly on the WGS84 spheroid. As an example I want to calculate the distance (in meters) between two linestrings. The postgresql code would be


SELECT
ST_Distance(
ST_GeomFromText('LINESTRING(13.45 52.47,13.46 52.48)')::Geography,
ST_GeomFromText('LINESTRING(13.00 52.00,13.1 52.2)')::Geography
)

The output is 38364.138314182


With geotools I would imagine something like



import org.geotools.geometry.jts.JTSFactoryFinder
import org.geotools.geometry.jts.WKTReader2,
import org.geotools.factory.Hints,
import org.geotools.referencing.crs.DefaultGeographicCRS

Hints hints = new Hints(Hints.CRS, DefaultGeographicCRS.WGS84)
GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(hints)
WKTReader2 reader = new WKTReader2(geometryFactory)
Geometry g1 = reader.read('LINESTRING(13.45 52.47,13.46 52.48)')
Geometry g2 = reader.read('LINESTRING(13.00 52.00,13.1 52.2)')

System.out.println(g1.distance(g2))

Now the output is 0.442040722106004 which is the same as


SELECT
ST_Distance(
ST_GeomFromText('LINESTRING(13.45 52.47,13.46 52.48)'),
ST_GeomFromText('LINESTRING(13.00 52.00,13.1 52.2)')
)

which is the euclidean distance in degrees. For points I can use the OrthodromicDistanceCalculator but not for lines. Plus, it is very inconvenient to specify that distance calculations should be based on WGS84 over and over again.




Answer



Using the following code I can get an answer of 38355.3256m (which is probably close enough) - I imagine that with a more careful choice of central point a more "accurate" result might be obtained.


    Geometry g1 = reader.read("LINESTRING(13.45 52.47,13.46 52.48)");
Geometry g2 = reader.read("LINESTRING(13.00 52.00,13.1 52.2)");
System.out.println(g1.distance(g2));
CoordinateReferenceSystem auto = CRS.decode("AUTO:42001,13.45,52.3");
System.out.println(auto);
MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84,
auto);
Geometry g3 = JTS.transform(g1, transform);

Geometry g4 = JTS.transform(g2, transform);
System.out.println(g3.distance(g4));

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...