How to Buffer a Point (Geography) by 1 Meter in Entity Framework Core with Net Topology Suite

Question

I'm using Entity Framework Core (3.1.0) with the Net Topology Suite packages as is recommended here, and I'm having a very hard time buffering a point to create a circle polygon around it. If I pass a value of 1 to Buffer(), where I'm under the belief that it means 1 meter, I end up with a circle with a radius close to 120 miles.

What value should I use to represent 1 meter? I'm using a GeometryFactory with an SRID of 4326 to create the points and polygons. Here's roughly the code I'm using:

GeometryFactory Geography = NtsGeometryServices.Instance.CreateGeometryFactory(4326);

var point = Geography.CreatePoint(coordinate.Latitude, coordinate.Longitude);
var polygonCoordinates = point.Buffer(1).Normalized().Reverse().Coordinates;
var polygon = Geography.CreatePolygon(polygonCoordinates);

Answer 1

After I posted my question, I decided to read through the Spatial docs for EF Core again. As @Michael Entin mentioned in the comments, the docs state that on the client side NTS ignores the SRID when it performs calculations. Your input has to be projected into a different coordinate system, calculated, and projected back. In my case I had to go from 4326 to 2855.

NOTE: 2855 is the US transformation in meters, for feet use 2926. For spatial references look here, and for their visualizations look here.

I looked through the example code, and decided to copy it into LINQPad so I can go through it in more detail. As it exists currently, the sample code "works" but only on a geometry with a single coordinate. If you pass in a polygon geometry which has multiple coordinates (such as the buffer polygon I'm using), then it just transforms the first one and ignores the rest. You have to loop through the remaining coordinates, make Points out of them, transform the points, and get the coordinate out of them.

After a bit more trial and error, I updated the example code so it can handle single and multi coordinate geometries. Checking the transformed coordinates on Google Maps seems to verify that it's correct.

Here's the updated projection code for anyone else struggling with this:

public static class CoordinateExtensions {
    public static Coordinate ProjectTo(
        this Coordinate coordinate,
        int fromSrid,
        int toSrid) {
        var point = new Point(coordinate) {
            SRID = fromSrid
        };

        return point.ProjectTo(toSrid).Coordinate;
    }
}

public static class GeometryExtensions {
    private static readonly CoordinateSystemServices Services = new CoordinateSystemServices(new Dictionary<int, string> {
        { 4326, GeographicCoordinateSystem.WGS84.WKT },
        { 2855, @"
            PROJCS[""NAD83(HARN) / Washington North"",
                GEOGCS[""NAD83(HARN)"",
                    DATUM[""NAD83_High_Accuracy_Reference_Network"",
                        SPHEROID[""GRS 1980"",6378137,298.257222101,
                            AUTHORITY[""EPSG"",""7019""]],
                        TOWGS84[0, 0, 0, 0, 0, 0, 0],
                        AUTHORITY[""EPSG"", ""6152""]],
                    PRIMEM[""Greenwich"", 0,
                        AUTHORITY[""EPSG"", ""8901""]],
                    UNIT[""degree"", 0.0174532925199433,
                        AUTHORITY[""EPSG"", ""9122""]],
                    AUTHORITY[""EPSG"", ""4152""]],
                PROJECTION[""Lambert_Conformal_Conic_2SP""],
                PARAMETER[""standard_parallel_1"", 48.73333333333333],
                PARAMETER[""standard_parallel_2"", 47.5],
                PARAMETER[""latitude_of_origin"", 47],
                PARAMETER[""central_meridian"", -120.8333333333333],
                PARAMETER[""false_easting"", 500000],
                PARAMETER[""false_northing"", 0],
                UNIT[""metre"", 1,
                    AUTHORITY[""EPSG"", ""9001""]],
                AXIS[""X"", EAST],
                AXIS[""Y"", NORTH],
                AUTHORITY[""EPSG"", ""2855""]]
        " }
    });

    public static Geometry ProjectTo(
        this Geometry geometry,
        int toSrid) {
        if (geometry.Coordinates.Length == 1) {
            return geometry.ProjectToSingle(toSrid);
        }

        return geometry.ProjectToMany(toSrid);
    }

    private static Geometry ProjectToSingle(
        this Geometry geometry,
        int toSrid) {
        var transformation = Services.CreateTransformation(geometry.SRID, toSrid);
        var transformer = new MathTransformFilter(transformation.MathTransform);
        var transformed = geometry.Copy();

        transformed.Apply(transformer);
        transformed.SRID = toSrid;

        return transformed;
    }

    private static Geometry ProjectToMany(
        this Geometry geometry,
        int toSrid) {
        var fromSrid = geometry.SRID;
        var transformed = geometry.Copy();

        for (var i = 0; i < transformed.Coordinates.Length; i++) {
            var coordinate = transformed.Coordinates[i].ProjectTo(fromSrid, toSrid);

            transformed.Coordinates[i].CoordinateValue = coordinate.CoordinateValue;
        }

        transformed.SRID = toSrid;

        return transformed;
    }
}

public sealed class MathTransformFilter :
    ICoordinateSequenceFilter {
    private readonly MathTransform Transform;

    public MathTransformFilter(
        MathTransform transform) => Transform = transform;

    public bool Done => false;
    public bool GeometryChanged => true;

    public void Filter(
        CoordinateSequence seq,
        int i) {
        var result = Transform.Transform(new[] {
            seq.GetOrdinate(i, Ordinate.X),
            seq.GetOrdinate(i, Ordinate.Y)
        });

        seq.SetOrdinate(i, Ordinate.X, result[0]);
        seq.SetOrdinate(i, Ordinate.Y, result[1]);
    }
}