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