1

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?

2
  • 1
    Web Mercator is notorious for distance inaccuracy. You should not be using 3857 to calculate buffers. Instead use a local appropriate UTM or Equal Area projection.
    – Vince
    Oct 24, 2020 at 18:47
  • Equal Area CRS codes (54012, 54009) seems to be not supported with Java's GeoTools. I am not sure if I am able to use local UTM, because my goal is to calculate such polygon between different zones. I wonder if it is possible to adjust my result to real distance somehow.
    – staszko032
    Oct 24, 2020 at 19:11

2 Answers 2

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.

6
  • This is exactly what I did initially and it works perfect for like European region. Thank you for that. The problem appears when I tried to buffer a line between Moscow and New York and it said org.geotools.referencing.operation.projection.ProjectionException: The transform result may be 4,780,472.998 meters away from the expected position. Are you sure that the input coordinates are inside this map projection area of validity? The point is located 45°57.4'W away from the central meridian and 70°31.0'N away from the latitude of origin. The projection is "Transverse_Mercator".
    – staszko032
    Oct 25, 2020 at 17:13
  • Then I tried EPSG:4326 -> EPSG:3857 transformation first, to at least get the result.
    – staszko032
    Oct 25, 2020 at 17:15
  • which is why I used the phrase - rough and ready. 3857 may give a "result" but it is very worng.
    – Ian Turton
    Oct 26, 2020 at 8:36
  • Then what do you think I should do? I need to calculate such buffers in meters within different cities in the world, I have lat long coordinates and then I need such polygons in MongoDB queries, so their format must be normalized degrees (-180/180, -90,90), eventually I need to display them on openstreetmap. I am out of ideas.
    – staszko032
    Oct 26, 2020 at 11:11
  • 1
    consider using PostGIS or port their Geography Buffer code to GeoTools
    – Ian Turton
    Oct 26, 2020 at 11:29
0

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;
    }

}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.