Transforming point using NetTopologySuite

Question

I'm exploring how we can use NetTopologySuite to address a geoprocessing task. based on some code sample I found, I was able to piece together a sample application that reads a shapefile and checks if a point falls within a polygon. The issue is, my point is in WGS84 but the shapefile is in EPSG:2810. I need to transform the input point to the CRS of the shapefile.

How do I convert a point from WGS to EPSG:2810 using NetTopologySuite 1.13.3?

Here's some of the code I have so far:

string shpFilename = @"C:\...\polygons.shp";
GeometryFactory factory = new GeometryFactory();
var coord = new Coordinate(-90.00, 45.00);
var point = factory.CreatePoint(coord);

I should also note that when I manually convert from WGS84 to the projected system and use those converted values, my code works as designed. This confirms this is a matter of transforming the coordinates before attempting to determine if they fall within a polygon.

Answer 1

The standard way of making point-to-point coordinate conversions in NTS is by using the ProjNet4GeoAPI library.

Here is a coordinate transformation example taken from the unit tests:

public void TestTransformListOfCoordinates()
{
    CoordinateSystemFactory csFact = new CoordinateSystemFactory();
    CoordinateTransformationFactory ctFact = new CoordinateTransformationFactory();

    ICoordinateSystem utm35ETRS = csFact.CreateFromWkt(
            "PROJCS[\"ETRS89 / ETRS-TM35\",GEOGCS[\"ETRS89\",DATUM[\"D_ETRS_1989\",SPHEROID[\"GRS_1980\",6378137,298.257222101]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",27],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]");

    IProjectedCoordinateSystem utm33 = ProjectedCoordinateSystem.WGS84_UTM(33, true);

    ICoordinateTransformation trans = ctFact.CreateFromCoordinateSystems(utm35ETRS, utm33);

    Coordinate[] points = new Coordinate[]
    {
        new Coordinate(290586.087, 6714000), new Coordinate(290586.392, 6713996.224),
        new Coordinate(290590.133, 6713973.772), new Coordinate(290594.111, 6713957.416),
        new Coordinate(290596.615, 6713943.567), new Coordinate(290596.701, 6713939.485)
    };

    Coordinate[] tpoints = trans.MathTransform.TransformList(points).ToArray();
    for (int i = 0; i < points.Length; i++)
        Assert.That(tpoints[i].Equals(trans.MathTransform.Transform(points[i])));
}

Answer 2

The answer works fine with points. But if you do this with polygons with holes it will return points out of order. Then when you turn those points back into a polygon it will complain about LineStrings not closing.

I ended up using the ProjNet library to build my transformation and then NetTopologySuite has a utility that handles polygons.

Looks something like this. You can get the WKT projection text from SpatialReference.org. You could easily expand this for all geometry types.

static GeoAPI.Geometries.IPolygon ProjectGeometry(GeoAPI.Geometries.IPolygon polygeo, string FromWKT, string ToWKT)
{
    var SourceCoordSystem = new CoordinateSystemFactory().CreateFromWkt(FromWKT);
    var TargetCoordSystem = new CoordinateSystemFactory().CreateFromWkt(ToWKT);

    var trans = new CoordinateTransformationFactory().CreateFromCoordinateSystems(SourceCoordSystem, TargetCoordSystem);

    var poly = NetTopologySuite.CoordinateSystems.Transformations.GeometryTransform.TransformPolygon(polygeo.Factory, polygeo, trans.MathTransform);
    return poly;
}

Update
The above code doesn't work with the new version of NetTopologySuite. Below is an updated version that works with NetTopologySuite 2.0.0.

static NetTopologySuite.Geometries.Geometry ProjectGeometry(NetTopologySuite.Geometries.Geometry geom, string FromWKT, string ToWKT)
{
    var SourceCoordSystem = new CoordinateSystemFactory().CreateFromWkt(FromWKT);
    var TargetCoordSystem = new CoordinateSystemFactory().CreateFromWkt(ToWKT);

    var trans = new CoordinateTransformationFactory().CreateFromCoordinateSystems(SourceCoordSystem, TargetCoordSystem);

    var projGeom = Transform(geom, trans.MathTransform);

    return projGeom;
}

static NetTopologySuite.Geometries.Geometry Transform(NetTopologySuite.Geometries.Geometry geom, MathTransform transform)
{
    geom = geom.Copy();
    geom.Apply(new MTF(transform));
    return geom;
}

sealed class MTF : NetTopologySuite.Geometries.ICoordinateSequenceFilter
{
    private readonly MathTransform _mathTransform;

    public MTF(MathTransform mathTransform) => _mathTransform = mathTransform;

    public bool Done => false;
    public bool GeometryChanged => true;
    public void Filter(NetTopologySuite.Geometries.CoordinateSequence seq, int i)
    {
        double x = seq.GetX(i);
        double y = seq.GetY(i);
        double z = seq.GetZ(i);
        _mathTransform.Transform(ref x, ref y, ref z);
        seq.SetX(i, x);
        seq.SetY(i, y);
        seq.SetZ(i, z);
    }
}

Answer 3

Faced the same need. I also found an example from the test. They revealed the details. But this approach is not convenient when working with the IGeometry type (IPoint, ILineString, etc).

After some research, I found how to transform geometry in a more convenient way (Thanks to NTS github repository :)

1) we should create transformation:

ICoordinateTransformation Wgs84ToWebMercator =
        (new CoordinateTransformationFactory())
        .CreateFromCoordinateSystems(GeographicCoordinateSystem.WGS84, 
                ProjectedCoordinateSystem.WebMercator);

2) declare filter, which will iterate through geometry sequence and made all work:

public class CoordinateTransformationFilter : ICoordinateSequenceFilter
{
    private readonly ICoordinateTransformation _transformation;

    public CoordinateTransformationFilter(ICoordinateTransformation transformation)
    {
        _transformation = transformation ??
            throw new ArgumentNullException($"{nameof(transformation)} can't be null.");
    }

    /// <summary>
    /// Yes, always true.
    /// </summary>
    public bool Done => true;

    /// <summary>
    /// Automatic call IGeometry.GeometryChanged() method after tranformation.
    /// </summary>
    public bool GeometryChanged => true;

    public void Filter(ICoordinateSequence seq, int i)
    {
        _transformation.MathTransform.Transform(seq);
    }
}

3) create filter:

var transformFilter = new CoordinateTransformationFilter(Wgs84ToWebMercator );

4) Profit:

// Get point from anywhere in WGS84 (aka EPSG:4326)
IPoint srcPoint = ...;
// You can skip it. But your geometry will be changed.
IPoint projectedPoint = srcPoint.Clone() as IPoint;
// Magic happen here. Our point in WebMercator (aka EPSG:3857)
projectedPoint.Apply(transformFilter);

UPDATE 2019-07-09: this solution for ProjNET4GeoAPI v1.4.1. I assume that it will not work with new version.