Buffer line between two coordinates with meters

Question

I am trying to build a polygon (with long/lat coordinates) over a line between two points (long/lat as well). Basically what I used is line buffering with meters as an distance.

However I have noticed that my polygon is not exactly created with choosen buffer offset (e.g 100km as buffer gives me like ~75km range polygon on OpenStreetMap). How can I achieve correct precision of buffer in meters?

I am using Java GeoTools.

Line builder:

public GeoLineBuilder() throws SchemaException {
    this.lineSchema = DataUtilities.createType(
            "",
            "Location",
            "locations:LineString:srid=4326," + "id:Integer"
    );
}

SimpleFeature buildLine(List<LocationPoint> points) {
    SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(lineSchema);
    GeometryFactory geometryFactory = new GeometryFactory();
    List<Coordinate> coordinates = points.stream().map(point -> new Coordinate(point.getLongtitude(), point.getLatitude())).collect(Collectors.toList());
    LineString line = geometryFactory.createLineString(coordinates.toArray(new Coordinate[] {}));
    featureBuilder.add(line);
    return featureBuilder.buildFeature(null);
}

Polygon over the line with buffer in meters:

private SimpleFeature bufferLine(SimpleFeature feature, long distance) throws Exception {
    GeometryAttribute gProp = feature.getDefaultGeometryProperty();
    CoordinateReferenceSystem origCRS = gProp.getDescriptor()
            .getCoordinateReferenceSystem();

    Geometry geom = (Geometry) feature.getDefaultGeometry();
    Geometry pGeom = geom;
    MathTransform toTransform, fromTransform = null;

    if (!(origCRS instanceof ProjectedCRS)) {

        Point c = geom.getCentroid();
        double x = c.getCoordinate().x;
        double y = c.getCoordinate().y;

        String code = "EPSG:3857";
        CoordinateReferenceSystem auto = CRS.decode(code, true);
        toTransform = CRS.findMathTransform(
                DefaultGeographicCRS.WGS84, auto);
        fromTransform = CRS.findMathTransform(auto,
                DefaultGeographicCRS.WGS84);
        pGeom = JTS.transform(geom, toTransform);
    }

    Geometry out = pGeom.buffer(distance);
    Geometry retGeom = out;

    if (!(origCRS instanceof ProjectedCRS)) {
            retGeom = JTS.transform(out, fromTransform);
    }

    SimpleFeatureType schema = feature.getFeatureType();
    SimpleFeatureTypeBuilder ftBuilder = new SimpleFeatureTypeBuilder();
    ftBuilder.setCRS(origCRS);
    ftBuilder.addAll(schema.getAttributeDescriptors());
    ftBuilder.setName(schema.getName());

    SimpleFeatureType nSchema = ftBuilder.buildFeatureType();
    SimpleFeatureBuilder builder = new SimpleFeatureBuilder(nSchema);
    List<Object> atts = feature.getAttributes();
    for (int i = 0; i < atts.size(); i++) {
        if (atts.get(i) instanceof Geometry) {
            atts.set(i, retGeom);
        }
    }
    return builder.buildFeature(null, atts.toArray());
}

Basically:

• I am building a line in EPSG:4326 to have a line with long/lat coordinates
• I am doing an EPSG:4326 -> EPSG:3857 transformation to be able to buffer by meters
• Actual buffering
• Reverse transformation to coordinates EPSG:4326 <- EPSG:3857
• Displaying received polygon on the openstreetmap map

I saw that when I choose Poznan and Paris coordinates as a line pair - received polygon was for sure smaller than selected buffering (100000 meters), it was like ~75000 meters instead.

How to approach this problem?

Answer 1

First DO NOT ever use Web Mercator for an application where you need distances to be correct.

You don't say what level of accuracy you need for your buffers or what the maximum lengths of line are so it's hard to give a precise answer. But for a relatively rough and ready solution I use the GeoTools AUTO projections (taken from the WMS specification) - see this StackExchange question and answer for a solution.

If you need more precision then you could look at using PostGIS and it's geography methods, Or even port those methods to GeoTools.

Answer 2

This is a different approach for creating a buffered line. It's gets thicker on the equator. Maybe you could work it around by adding intermediate points

import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.geotools.util.factory.Hints;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.MultiPolygon;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.util.GeometricShapeFactory;

public class Scratch {

    private static final GeometryFactory GEOMETRY_FACTORY = JTSFactoryFinder.getGeometryFactory(
            new Hints(Hints.CRS, DefaultGeographicCRS.WGS84));
    private static final int CIRCLE_SIDES = 16;
    private static final double LATITUDE_METERS_IN_ONE_DEGREE = 111_320;
    private static final double EQUATORIAL_CIRCUMFERENCE_OF_EARTH_IN_METERS = 40_075_000;
    private static final int DEGREES_FULL_CIRCLE = 360;

    public static void main(String[] args) {
        toMultiPolygon(Arrays.asList(
                toCircle(new Coordinate(52.2297, 21.0122), 1), //Warsaw
                toCircle(new Coordinate(-34.6037, -58.3816), 2) // Buenos Aires
        )).convexHull();

    }

    public static Polygon toCircle(Coordinate center, Number radiusInMeters) {
        GeometricShapeFactory shapeFactory = new GeometricShapeFactory();
        shapeFactory.setNumPoints(CIRCLE_SIDES);
        shapeFactory.setCentre(center);
        shapeFactory.setHeight(toHeightInRadians(radiusInMeters.doubleValue() * 2));
        shapeFactory.setWidth(toWidthInRadians(radiusInMeters.doubleValue() * 2, center.y));
        return shapeFactory.createEllipse();
    }

    public static MultiPolygon toMultiPolygon(List<Polygon> polygons) {
        return GEOMETRY_FACTORY.createMultiPolygon(
                polygons.toArray(Polygon[]::new));
    }

    private static double toWidthInRadians(Number size, double latitude) {
        return size.doubleValue() / (EQUATORIAL_CIRCUMFERENCE_OF_EARTH_IN_METERS * Math.cos(Math.toRadians(latitude)) / DEGREES_FULL_CIRCLE);
    }

    private static double toHeightInRadians(Number size) {
        return size.doubleValue() / LATITUDE_METERS_IN_ONE_DEGREE;
    }

}