Merging polygons on curved surface

Question

I tried the clipper library for merging polygons (c#) that I want to import into Google Earth. It worked with a flat-surface sample of two circles, but there are two problems with using longitute and latitude as X and Y:

1. Clipper only supports long as datatype (you can get around that be using a multiplier before passing the values to Clipper - there is already a question about that, but it doesn't address the second point
2. Clipper assumes the surface is flat, so the output coordinates are not where they are expected. It's not possible to just think of X and Y as longitude and latitude and add value x to them. It may work at the equator, but higher up north adding the same value to longitude would result in a shorter distance.

I was wondering if there are already implementations / libraries for that scenario.

Update: As the source code that was used has been requested, I created a small sample application that contains the relevant code:

internal class Polygon : List<DoublePoint> { }
internal class Polygons : List<List<DoublePoint>> { }

internal class IntPolygon : List<IntPoint> { }
internal class IntPolygons : List<List<IntPoint>> { }

class Program
{
    static void Main(string[] args)
    {
        MercatorProjection mercatorProjection = new MercatorProjection();

        List<GeoCoordinate> circle1Coordinates = GenerateCirclePolygonByMeters(new GeoCoordinate { Longitude = 50, Latitude = 50 }, 5000);
        List<GeoCoordinate> circle2Coordinates = GenerateCirclePolygonByMeters(new GeoCoordinate { Longitude = 50.05, Latitude = 50 }, 5000);

        Polygons polygons = new Polygons();

        Polygon circle1 = new Polygon();
        foreach(var coordinate in circle1Coordinates)
        {
            DoublePoint point = mercatorProjection.FromLatLngToPoint(coordinate.Latitude, coordinate.Longitude, 15);
            circle1.Add(point);
        }

        Polygon circle2 = new Polygon();
        foreach (var coordinate in circle2Coordinates)
        {
            DoublePoint point = mercatorProjection.FromLatLngToPoint(coordinate.Latitude, coordinate.Longitude, 15);
            circle2.Add(point);
        }

        polygons.Add(circle1);
        polygons.Add(circle2);

        // Workaround: get int-values
        IntPolygons intPolygons = ConvertToIntPolygons(polygons);

        // Merge
        IntPolygons mergedIntPolygons = MergePolygons(intPolygons);

        // Workaroud: Get doubles again
        Polygons mergedPolygons = ConvertToDoublePolygons(mergedIntPolygons);

        // Convert back to spherical surface
        // GeoCoordinate class is from System.Device
        List<GeoCoordinate> mergedCoordinates = new List<GeoCoordinate>();
        foreach (var polygon in mergedPolygons)
        {
            foreach (var point in polygon)
            {
                GeoCoordinate coordinate = mercatorProjection.FromPointToLatLng(point, 15);
                mergedCoordinates.Add(coordinate);
            }
        }

        // Generate output csv-list
        WriteOutputFile(mergedCoordinates);
    }

    private static void WriteOutputFile(List<GeoCoordinate> coordinates, string filename = "")
    {
        string uniquename = DateTime.Now.ToString("yyyyMMdd_HHmmss") + ".txt";
        string filenameToUse = String.IsNullOrEmpty(filename) ? uniquename : filename;
        var fileStream = File.OpenWrite(filenameToUse);

        using (var writer = new StreamWriter(fileStream))
        {
            // header
            writer.WriteLine("lat,long");

            // content
            foreach (var coordinate in coordinates)
            {
                writer.WriteLine(coordinate.Latitude.ToString(CultureInfo.InvariantCulture) + "," + coordinate.Longitude.ToString(CultureInfo.InvariantCulture));
            }

            writer.Close();
        }
    }

    private static Polygons ConvertToDoublePolygons(IntPolygons polygons)
    {
        double multiplier = 100;

        Polygons doublePolygons = new Polygons();
        foreach (var intPolygon in polygons)
        {
            Polygon doublePolygon = new Polygon();
            foreach (var intPoint in intPolygon)
            {
                doublePolygon.Add(new DoublePoint { X = intPoint.X / multiplier, Y = intPoint.Y / multiplier });
            }

            doublePolygons.Add(doublePolygon);
        }

        return doublePolygons;
    }

    private static IntPolygons ConvertToIntPolygons(Polygons polygons)
    {
        double multiplier = 100;
        IntPolygons intPolygons = new IntPolygons();
        foreach (var doublePolygon in polygons)
        {
            IntPolygon intPolygon = new IntPolygon();
            foreach (var doublePoint in doublePolygon)
            {
                intPolygon.Add(new IntPoint { X = (int)(doublePoint.X * multiplier), Y = (int)(doublePoint.Y * multiplier) });
            }

            intPolygons.Add(intPolygon);
        }

        return intPolygons;
    }

    private static List<GeoCoordinate> GenerateCirclePolygonByMeters(GeoCoordinate position, double distanceMeters)
    {
        // calc distance in coordinates-system
        double latitude2 = CalculateDerivedPosition(position, distanceMeters, 0).Latitude;
        double yRadius = latitude2 - position.Latitude;
        double longitude2 = CalculateDerivedPosition(position, distanceMeters, 90).Longitude;
        double xRadius = longitude2 - position.Longitude;

        var circle = GenerateCirclePolygon(position.Longitude, position.Latitude, xRadius, yRadius, 16);

        List<GeoCoordinate> circleCoordinates = new List<GeoCoordinate>();
        foreach(var point in circle)
        {
            circleCoordinates.Add(new GeoCoordinate { Latitude = point.Y, Longitude = point.X });
        }

        return circleCoordinates;
    }

    private static Polygon GenerateCirclePolygon(double xStart, double yStart, double xRadius, double yRadius, int points)
    {
        Polygon polygon = new Polygon();
        double slice = 2 * Math.PI / points;


        for (int i = 0; i < points; i++)
        {
            double angle = slice * i;
            Console.WriteLine(angle);
            double x = xRadius * Math.Cos(angle);
            double y = yRadius * Math.Sin(angle);

            polygon.Add(new DoublePoint(xStart + x, yStart + y));
        }

        return polygon;
    }

    // Source: https://stackoverflow.com/questions/1125144/how-do-i-find-the-lat-long-that-is-x-km-north-of-a-given-lat-long
    public static GeoCoordinate CalculateDerivedPosition(GeoCoordinate source, double rangeMeters, double bearing)
    {
        double radiansToDegrees = 57.2957795;
        double degreesToRadians = 0.0174532925;
        double twoPi = Math.PI * 2;
        double earthRadius = 6378137.0;

        double latA = source.Latitude * degreesToRadians;
        double lonA = source.Longitude * degreesToRadians;
        double angularDistance = rangeMeters / earthRadius;
        double trueCourse = bearing * degreesToRadians;

        double lat = Math.Asin(
            Math.Sin(latA) * Math.Cos(angularDistance) +
            Math.Cos(latA) * Math.Sin(angularDistance) * Math.Cos(trueCourse));

        double dlon = Math.Atan2(
            Math.Sin(trueCourse) * Math.Sin(angularDistance) * Math.Cos(latA),
            Math.Cos(angularDistance) - Math.Sin(latA) * Math.Sin(lat));

        double lon = ((lonA + dlon + Math.PI) % twoPi) - Math.PI;

        return new GeoCoordinate(
            lat * radiansToDegrees,
            lon * radiansToDegrees);
    }

    private static IntPolygons MergePolygons(IntPolygons polygons)
    {
        Clipper clipper = new Clipper();

        clipper.AddPaths(polygons, PolyType.ptSubject, true);

        IntPolygons mergedPolygons = new IntPolygons();
        clipper.Execute(ClipType.ctUnion, mergedPolygons,
            PolyFillType.pftNonZero, PolyFillType.pftNonZero);

        return mergedPolygons;
    }

Update 2: There was actually an error in the way that the circles were created. With the correction the circles get created and merged correctly:

I found out that I can use the Mercator Projection to get the points of the circles in the image above projected to a flat surface. Then merging could be done using Clipper. Finally all points were projected to the sphere again (the reverse of the Mercator Projection).

This would be a workaround (will be testing it in more depth). So the question is kind of answered, but it would be interesting if there other solutions with less workarounds.