Dalton Filho

The first and simplest algorithm I will mention deals with segmenting LineStrings so that each segment will have a fixed arbitrary length. Your task is to divide Montreal’s Grand Prix track in segments of 100 meters.

To begin working on this algorithm, you need data, which you can get from the KML (EPSG:4326) of the track you see below:

Maven dependencies and repositories

<properties>
    <project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
    <geotools.version>12-SNAPSHOT</geotools.version>
</properties>
<dependencies>
    <dependency>
        <groupId>commons-lang</groupId>
        <artifactId>commons-lang</artifactId>
        <version>2.6</version>
    </dependency>
    <dependency>
        <groupId>commons-io</groupId>
        <artifactId>commons-io</artifactId>
        <version>2.4</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-xml</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools.xsd</groupId>
        <artifactId>gt-xsd-kml</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-referencing</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-render</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-epsg-wkt</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>java3d</groupId>
        <artifactId>vecmath</artifactId>
        <version>1.3.1</version>
    </dependency>
</dependencies>

<repositories>
    <repository>
        <id>maven2-repository.dev.java.net</id>
        <name>Java.net repository</name>
        <url>http://download.java.net/maven/2 </url>
    </repository>
    <repository>
        <id>osgeo</id>
        <name>Open Source Geospatial Foundation Repository</name>
        <url>http://download.osgeo.org/webdav/geotools </url>
    </repository>
    <repository>
        <snapshots>
            <enabled>true</enabled>
        </snapshots>
        <id>opengeo</id>
        <name>OpenGeo Maven Repository</name>
        <url>http://repo.opengeo.org </url>
    </repository>
</repositories>

Importing the example KML into your application

	

import java.io.BufferedInputStream;
import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStream;
import java.net.URISyntaxException;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
import java.util.Iterator;
import java.util.LinkedList;
import java.util.List;

import javax.xml.parsers.ParserConfigurationException;

import org.geotools.kml.KMLConfiguration;
import org.geotools.referencing.CRS;
import org.geotools.referencing.GeodeticCalculator;
import org.geotools.xml.Parser;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.NoSuchAuthorityCodeException;
import org.opengis.referencing.operation.TransformException;
import org.xml.sax.SAXException;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.LineString;
import com.vividsolutions.jts.geom.Point;
import com.vividsolutions.jts.geom.PrecisionModel;

@SuppressWarnings("unchecked")
private static LineString readGeometryFromKML(String kmlPath) throws IOException, SAXException, 
    ParserConfigurationException {

    Parser parser = new Parser(new KMLConfiguration());
    InputStream stream = new BufferedInputStream(new FileInputStream(kmlPath));
    LineString track = null;
    
    try {
        SimpleFeature feature = (SimpleFeature) parser.parse(stream);
        Collection placemarks = (Collection) feature.getAttribute("Feature");
        track = (LineString) placemarks.iterator().next().getDefaultGeometry();
    } finally {
        stream.close();
    }
    
    return track;
}

Now you’re in possession of a LineString object which you will have to divide in segments of 100m. Predictably, there is no method in GeoTools to provide such a specific functionality. What that means is that you’re in for a lot of fun!

The two determinant steps to create this algorithm are:

1. Determining the distance from point A to B
2. Interpolating a point within the line so that you can segment the line at the precise point where the specified distance is reached

Step #1: determining distances in GeoTools

GeoTools has a handy class called GeodeticCalculator for calculating the orthodromic distance between two coordinates in meters. You could also use the orthodromicDistance() method of the JTS class. The most important thing to remember when working with GeodeticCalculator or JTS is the fact that JTS calculates the distance atomically, in just one method. GeodeticCalculator, on the other hand, will require 3 method calls to calculate the same distance. This has implications in multi-threading scenarios depending on how you plan to share the GeodeticCalculator object.

Using GeodeticCalculator

	

LineString track = readGeometryFromKML("MontrealTrack.kml");
GeodeticCalculator calculator = new GeodeticCalculator(CRS.decode("EPSG:4326")); // KML uses WGS84
Coordinate[] coordinates = track.getCoordinates();

for (int c = 0; c < coordinates.length - 1; c++) {
    Coordinate c1 = coordinates[c];
    Coordinate c2 = coordinates[c + 1];
    
    calculator.setStartingGeographicPoint(c1.x, c1.y);
    calculator.setDestinationGeographicPoint(c2.x, c2.y);
    
    System.out.printf("Segment [%s-%s]. Length: %sm %n", c, c + 1, calculator.getOrthodromicDistance());
}

Using JTS

	

LineString track = readGeometryFromKML("MontrealTrack.kml");
CoordinateReferenceSystem crsWGS84 = CRS.decode("EPSG:4326");
Coordinate[] coordinates = track.getCoordinates();

for (int c = 0; c < coordinates.length - 1; c++) {
    Coordinate c1 = coordinates[c];
    Coordinate c2 = coordinates[c + 1];
    
    System.out.printf("Segment [%s-%s]. Length: %sm %n", c, c + 1, JTS.orthodromicDistance(c1, c2, crsWGS84));
}

Step #2: interpolating a point within a line

If you can calculate the length of a segment, you can easily sum up the lengths of consecutive segments until you reach a length that is equal or above the specified length — 100m in our case. Your algorithm will eventually reach a segment where its length + the accumulated length is above the required length. You now have to divide your segment at the precise point where you reach 100m. Once again, GeoTools won’t give you this specific functionality out of the box, so you will have to interpolate this point on your own.

One way to know this point is by using the properties of similar triangles.

Let a-b be your line segment. You can easily obtain its 3 important properties:

1. The distance a-b in the X axis (b.x – a.x). Let’s call it “DX”
2. The distance a-b in the Y axis (b.y – a.y). Let’s call it “DY”
3. The distance a-b in meters (Using GeodeticCalculator or JTS). Let’s call it “D”

Once in possession of these 3 properties you can form a triangle and treat your line segment as your hypotenuse.

When you divide a line, you can form another triangle in the same fashion, thus creating a child triangle with the same proportions.

Now you have everything you need to find the longitude (x) and latitude (y) of the point of segmentation (b):

The distance (d) is given by subtracting the required length from the accumulated length. If the required length is 100m and the accumulated length so far is 115m, d = 15m.

Implementation

The code below implements the steps discussed above. Adapt it and optimize it as you see fit!

List<LineString> createSegments(Geometry track, double segmentLength) throws NoSuchAuthorityCodeException, 
    FactoryException {

    GeodeticCalculator calculator = new GeodeticCalculator(CRS.decode("EPSG:4326")); // KML uses WGS84
    GeometryFactory geometryFactory = new GeometryFactory(new PrecisionModel(PrecisionModel.FLOATING), 4326);
    
    LinkedList<Coordinate> coordinates = new LinkedList<Coordinate>();
    Collections.addAll(coordinates, track.getCoordinates());
    
    double accumulatedLength = 0;
    List<Coordinate> lastSegment = new ArrayList<Coordinate>();
    List<LineString> segments = new ArrayList<LineString>();
    Iterator<Coordinate> itCoordinates = coordinates.iterator();
    
    for (int i = 0; itCoordinates.hasNext() && i < coordinates.size() - 1; i++) {
        Coordinate c1 = coordinates.get(i);
        Coordinate c2 = coordinates.get(i + 1);
        
        lastSegment.add(c1);
        
        calculator.setStartingGeographicPoint(c1.x, c1.y);
        calculator.setDestinationGeographicPoint(c2.x, c2.y);
        
        double length = calculator.getOrthodromicDistance();
        
        if (length + accumulatedLength >= segmentLength) {
            double offsetLength = segmentLength - accumulatedLength;
            double ratio = offsetLength / length;
            double dx = c2.x - c1.x;
            double dy = c2.y - c1.y;
            
            Coordinate segmentationPoint = new Coordinate(c1.x + (dx * ratio), 
                                                          c1.y + (dy * ratio)); 
            
            lastSegment.add(segmentationPoint); // Last point of the segment is the segmentation point
            segments.add(geometryFactory.createLineString(lastSegment.toArray(new Coordinate[lastSegment.size()])));
            
            lastSegment = new ArrayList<Coordinate>(); // Resets the variable since a new segment will be built
            accumulatedLength = 0D;
            coordinates.add(i + 1, segmentationPoint); 
        } else {
            accumulatedLength += length;
        }
    }
    
    lastSegment.add(coordinates.getLast()); // Because the last one is never added in the loop above
    segments.add(geometryFactory.createLineString(lastSegment.toArray(new Coordinate[lastSegment.size()])));
    
    return segments;
}

Result

Below you can see the end result of segmenting the Montreal Grand Prix track in segments of 100m. Note how the last segment in blue, measuring 26m, is the only segment measuring less than 100m since the total track length is not a multiple of 100. You can download this kml.