System.ArgumentException: Transformation cannot be computed at the poles.

[presto41]

I'm trying to transform coordinates using ProjNET 2.0.0 as nuget package in my project. Unfortunately I've got an exception during this operation.

The exception is throwing because of line 131:

/* Forward equations */
if (Math.Abs(Math.Abs(dLatitude) - HALF_PI) <= EPSLN)
    throw new ArgumentException("Transformation cannot be computed at the poles.");

in method RadiansToMeters(ref double lon, ref double lat) from ProjNet\CoordinateSystems\Projections\Mercator.cs

The source of coordinate system is EPSG:3879 and the target is EPSG:3857. The most interesting is fact that everything is working well when I'm using in my project implementation with ProjNET4GeoAPI v1.4.1 from nuget packages.

I prepared a test case code (using ProjNET 2.0.0):

static void Main(string[] args)
{
    ICoordinateTransformation transform;

    var csfact = new CoordinateSystemFactory(); 

    var wkt3857 = "PROJCS[\"WGS 84 / Pseudo - Mercator\",GEOGCS[\"Popular Visualisation CRS\",DATUM[\"Popular_Visualisation_Datum\",SPHEROID[\"Popular Visualisation Sphere\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7059\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6055\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4055\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Mercator_1SP\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"semi_minor\",6378137],PARAMETER[\"central_meridian\",0],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],AUTHORITY[\"EPSG\",\"3785\"],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH]]";
    var wkt3879 = "PROJCS[\"ETRS89 / ETRS - GK25FIN\",GEOGCS[\"ETRS89\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6258\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4258\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",25],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",25500000],PARAMETER[\"false_northing\",0],AUTHORITY[\"EPSG\",\"3879\"],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]";

    var crdSys3857 = csfact.CreateFromWkt(wkt3857);
    var crdSys3879 = csfact.CreateFromWkt(wkt3879);

    var fact = new CoordinateTransformationFactory();
    transform = fact.CreateFromCoordinateSystems(crdSys3879, crdSys3857);

    double x = 262571166.69437805;
    double y = 336848165.97773886;
    double z = Double.NaN;

    try
    {
        transform.MathTransform.Transform(ref x, ref y, ref z);
    }
    catch (Exception ex)
    {
        Console.WriteLine(ex);
    }
    Console.ReadLine();
}

Would you help me with that? It's seems to be a bug in ProjNET.

[airbreather]

The most interesting is fact that everything is working well when I'm using in my project implementation with ProjNET4GeoAPI v1.4.1 from nuget packages.

This test case throws the same exception in 1.4.1:

ConsoleApp0.csproj

<Project Sdk="Microsoft.NET.Sdk">

  <PropertyGroup>
    <OutputType>Exe</OutputType>
    <TargetFramework>netcoreapp3.1</TargetFramework>
  </PropertyGroup>

  <ItemGroup>
    <PackageReference Include="ProjNET4GeoAPI" Version="1.4.1" />
  </ItemGroup>

</Project>

Program.cs

using System;

using GeoAPI.CoordinateSystems.Transformations;

using ProjNet.CoordinateSystems;
using ProjNet.CoordinateSystems.Transformations;

static class Program
{
    static void Main(string[] args)
    {
        ICoordinateTransformation transform;

        var csfact = new CoordinateSystemFactory();

        var wkt3857 = "PROJCS[\"WGS 84 / Pseudo - Mercator\",GEOGCS[\"Popular Visualisation CRS\",DATUM[\"Popular_Visualisation_Datum\",SPHEROID[\"Popular Visualisation Sphere\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7059\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6055\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4055\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Mercator_1SP\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"semi_minor\",6378137],PARAMETER[\"central_meridian\",0],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],AUTHORITY[\"EPSG\",\"3785\"],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH]]";
        var wkt3879 = "PROJCS[\"ETRS89 / ETRS - GK25FIN\",GEOGCS[\"ETRS89\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6258\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4258\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",25],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",25500000],PARAMETER[\"false_northing\",0],AUTHORITY[\"EPSG\",\"3879\"],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]";

        var crdSys3857 = csfact.CreateFromWkt(wkt3857);
        var crdSys3879 = csfact.CreateFromWkt(wkt3879);

        var fact = new CoordinateTransformationFactory();
        transform = fact.CreateFromCoordinateSystems(crdSys3879, crdSys3857);

        double x = 262571166.69437805;
        double y = 336848165.97773886;
        double z = Double.NaN;

        try
        {
            transform.MathTransform.Transform(new[] { x, y, z });
        }
        catch (Exception ex)
        {
            Console.WriteLine(ex);
        }
        Console.ReadLine();
    }
}

[presto41]

Hmm, so I'm not sure why it doesn't work well. My implementation of TransformGeometry (where I'm using 1.4.1 version) looks like that:

public IGeometry TransformGeometry(IGeometry geometry)
{
    if (geometry == null)
    {
        return null;
    }
    var udata = geometry.UserData;
    var rgeom = default(IGeometry);
    if (geometry is IGeometryCollection geometryColl)
    {
        rgeom = GeometryTransform.TransformGeometryCollection(this.gfact, geometryColl, this.transform.MathTransform);
    }
    else
    {
        rgeom = GeometryTransform.TransformGeometry(this.gfact, geometry, this.transform.MathTransform);
    }
    rgeom.UserData = udata;
    return rgeom;
}

and it has been changed when I updated my ProjNET version and now the implementation looks like that:

public Geometry TransformGeometry(Geometry geometry)
{
    if (geometry == null)
    {
        return null;
    }
    var udata = geometry.UserData;
    var rgeom = Transform(geometry, this.transform.MathTransform);
    rgeom.UserData = udata;
    return rgeom;
}

public static T Transform<T>(T geom, MathTransform transform) where T : Geometry
{
    geom = (T)geom.Copy();
    geom.Apply(new MTF(transform));
    return geom;
}

private sealed class MTF : ICoordinateSequenceFilter
{
    private readonly MathTransform mathTransform;

    public MTF(MathTransform mathTransform) => this.mathTransform = mathTransform;

    public bool Done => false;
    public bool GeometryChanged => true;
    public void Filter(CoordinateSequence seq, int i)
    {
        var (x, y, z) = this.mathTransform.Transform(seq.GetX(i), seq.GetY(i), seq.GetZ(i));
        seq.SetX(i, x);
        seq.SetY(i, y);
        seq.SetZ(i, z);

    }
}

This changes are exactly as you were suggesting here

The geometry which I'm trying to transform is this same in both cases. So I really don't know where is the problem if ProjNET is working well.

mathTransform looks like in my test case.

[presto41]

I know now what the problem is. In old way (ProjNET4GeoAPI v1.4.1) when I wanted to transform a polygon I had to use a TransformPolygon(...) method from NetTopologySuite.CoordinateSystems.Transformations. If this method couldn't to create a polygon then was returning empty polygon. In new version of ProjNET v2.0.0 in this same situation is throwing an exception instead of return empty polygon. It was the reason why my code didn't work well. I didn't expect a exception here.