garykindel / ShapefileProjectionDemo

ShapefileProjectionDemo
1 stars 1 forks source link

What's next? #2

Closed pauldendulk closed 5 years ago

pauldendulk commented 5 years ago

I was able to start the application and load the states of the US. It shows correctly. So the custom projection is working properly in that case. Is that correct?

For there other cases it does not work correctly? What are you trying to achieve? Would you like a general solution that works for all projection files that are loaded on the fly?

garykindel commented 5 years ago

The States file uses WGS84 so it projects fine. The problem is the NY townships (NYTOWNS_POLY.shp) test case. It is not working properly when projected with Open Street maps. The test app has code to load the NY townships with and without Open street maps.

I would like a solution that can specifically load NY Townships properly using its prj file with Open Street maps ( using EPSG:3857). If I had this code, I'm sure I could use it as a pattern to handle other 'custom' projections.

For context, I have an existing geology/mining database app written in C#, winforms with SharpMap. I can display NY townships with Open Streets map using SharpMap in this app. Admittedly, I find projections, GeoAPI and ProjNet all abit confusing. The Projectionsample.cs in Mapsui helped me get this far. I'm rewriting this app in C#, WPF, MVVM with MapSui. I have a large amount of mining and geology data in shapefiles I would like to use with my WPF app.

BTW, I appreciate you taking the time to look into this projection issue. MapSui is an excellent open source project and it is performing well (as part of my WPF app) on my i7 windows 10 laptop with SSD.

pauldendulk commented 5 years ago

I made some changes to the CustomMinimalTransformation

using GeoAPI.CoordinateSystems;
using GeoAPI.CoordinateSystems.Transformations;
using Mapsui.Geometries;
using Mapsui.Projection;
using ProjNet.CoordinateSystems.Transformations;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using ProjNet.CoordinateSystems;

namespace WpfTestApp
{
    public class CustomMinimalTransformation : ITransformation
    {
        private readonly IDictionary<string, Func<double, double, Point>> _toLonLat = new Dictionary<string, Func<double, double, Point>>();
        private readonly IDictionary<string, Func<double, double, Point>> _fromLonLat = new Dictionary<string, Func<double, double, Point>>();

        public CustomMinimalTransformation()
        {
            _toLonLat["EPSG:4326"] = (x, y) => new Point(x, y);
            _fromLonLat["EPSG:4326"] = (x, y) => new Point(x, y);
            _toLonLat["EPSG:3857"] = SphericalMercator.ToLonLat;
            _fromLonLat["EPSG:3857"] = SphericalMercator.FromLonLat;
            _toLonLat["EPSG:CUSTOM"] = CustomProjectionToLonLat;
            _fromLonLat["EPSG:CUSTOM"] = CustomProjectionFromLonLat;
        }

        public IGeometry Transform(string fromCRS, string toCRS, IGeometry geometry)
        {
            Transform(geometry.AllVertices(), _toLonLat[fromCRS]);
            Transform(geometry.AllVertices(), _fromLonLat[toCRS]);
            return geometry; // this method should not have a return value
        }

        private static void Transform(IEnumerable<Point> points, Func<double, double, Point> transformFunc)
        {
            foreach (var point in points)
            {
                var transformed = transformFunc(point.X, point.Y);
                point.X = transformed.X;
                point.Y = transformed.Y;
            }
        }

        public BoundingBox Transform(string fromCRS, string toCRS, BoundingBox boundingBox)
        {
            Transform(boundingBox.AllVertices(), _toLonLat[fromCRS]);
            Transform(boundingBox.AllVertices(), _fromLonLat[toCRS]);
            return boundingBox; // this method not have a return value
        }

        public bool? IsProjectionSupported(string fromCRS, string toCRS)
        {
            if (!_toLonLat.ContainsKey(fromCRS)) return false;
            if (!_fromLonLat.ContainsKey(toCRS)) return false;
            return true;
        }

        CoordinateTransformationFactory _ctFac;
        ICoordinateTransformation _ctTo;   //Custom to WGS84
        ICoordinateTransformation _ctFrom; //WGS84 to Custom

        public void LoadSourceWKT(string filepath)
        {
            //@"C:\DRC_Data\Arcview\USA\Townships\NYTOWNS_POLY.prj";

            ICoordinateSystemFactory csFac = new ProjNet.CoordinateSystems.CoordinateSystemFactory();
            string file = @"NYTOWNS_POLY.prj";
            string wkt = System.IO.File.ReadAllText(file);
            ICoordinateSystem csFrom = csFac.CreateFromWkt(wkt);
            _ctFac = new CoordinateTransformationFactory();

            _ctTo = _ctFac.CreateFromCoordinateSystems(
                csFrom,
                GeographicCoordinateSystem.WGS84);

            _ctFrom = _ctFac.CreateFromCoordinateSystems(
                GeographicCoordinateSystem.WGS84, 
                csFrom);
        }

        public Point CustomProjectionFromLonLat(double lon, double lat)
        {
            Point pt = new Point(lon, lat);                       
            double[] result= _ctFrom.MathTransform.Transform(pt.ToDoubleArray());
            return new Point(result[0], result[1]);
        }

        public Point CustomProjectionToLonLat(double x, double y)
        {
            Point pt = new Point(x, y);
            double[] result = _ctTo.MathTransform.Transform(pt.ToDoubleArray());
            return new Point(result[0], result[1]);
        }

    }
}
pauldendulk commented 5 years ago

Performance is not great for shapefiles though. I usually use the RasterizingLayer but this does not work for transformations. The workaround would be to create a custom provider which has the shapefile provider as member an transforms the coordinates before returning them to the layer. Look at the Mapsui implementation of transformation on how to do this.

garykindel commented 5 years ago

May projection is / was wrong. I replaced calls to SphericalMercator with transform in CustomMinimalTransformation. It initially projected lower 48 states to below southern edge of western Africa. In CustomMinimalTransformation I changed the following: `public class CustomMinimalTransformation : ITransformation { private readonly IDictionary<string, Func<double, double, Point>> _toLonLat = new Dictionary<string, Func<double, double, Point>>(); private readonly IDictionary<string, Func<double, double, Point>> _fromLonLat = new Dictionary<string, Func<double, double, Point>>();

    public CustomMinimalTransformation()
    {
        _toLonLat["EPSG:4326"] = (x, y) => new Point(x, y);
        _fromLonLat["EPSG:4326"] = (x, y) => new Point(x, y);
         _toLonLat["EPSG:3857M"] = SphericalMercator.ToLonLat;
         _fromLonLat["EPSG:3857M"] = SphericalMercator.FromLonLat;
         _toLonLat["EPSG:3857"] = CustomProjectionToLonLat;
         _fromLonLat["EPSG:3857"] = CustomProjectionFromLonLat;
        _toLonLat["EPSG:CUSTOM"] = CustomProjectionToLonLat;
        _fromLonLat["EPSG:CUSTOM"] = CustomProjectionFromLonLat;
    }

    public IGeometry Transform(string fromCRS, string toCRS, IGeometry geometry)
    {
        Transform(geometry.AllVertices(), _toLonLat[fromCRS]);
        Transform(geometry.AllVertices(), _fromLonLat[toCRS]);
        return geometry; // this method should not have a return value
    }

    private static void Transform(IEnumerable<Point> points, Func<double, double, Point> transformFunc)
    {
        foreach (var point in points)
        {
            var transformed = transformFunc(point.X, point.Y);
            point.X = transformed.X;
            point.Y = transformed.Y;
        }
    }

    public BoundingBox Transform(string fromCRS, string toCRS, BoundingBox boundingBox)
    {
        Transform(boundingBox.AllVertices(), _toLonLat[fromCRS]);
        Transform(boundingBox.AllVertices(), _fromLonLat[toCRS]);
        return boundingBox; // this method not have a return value
    }

    public bool? IsProjectionSupported(string fromCRS, string toCRS)
    {
        if (!_toLonLat.ContainsKey(fromCRS)) return false;
        if (!_fromLonLat.ContainsKey(toCRS)) return false;
        return true;
    }

    CoordinateTransformationFactory _ctFac;
    ICoordinateTransformation _ctTo;   
    ICoordinateTransformation _ctFrom; 

    public void LoadSourceWKT(string filepath)
    {
        ICoordinateSystemFactory csFac = new ProjNet.CoordinateSystems.CoordinateSystemFactory();
        ICoordinateSystem csSource = null;
        ICoordinateSystem csTarget = null;
        if (!String.IsNullOrWhiteSpace(filepath) && System.IO.File.Exists(filepath))
        {
            string wkt = System.IO.File.ReadAllText(filepath);
            csSource = csFac.CreateFromWkt(wkt);
            csTarget = ProjectedCoordinateSystem.WebMercator;
        }
        else
        {                
            csSource = GeographicCoordinateSystem.WGS84;
            csTarget = ProjectedCoordinateSystem.WebMercator;
        }
        _ctFac = new CoordinateTransformationFactory();

        _ctTo = _ctFac.CreateFromCoordinateSystems(csSource, csTarget);
        _ctFrom = _ctFac.CreateFromCoordinateSystems(csTarget, csSource);
    }

    public Point CustomProjectionFromLonLat(double lon, double lat)
    {
        Point pt = new Point(lon, lat);                       
        double[] result= _ctTo.MathTransform.Transform(pt.ToDoubleArray());
        return new Point(result[0], result[1]);
    }

    public Point CustomProjectionToLonLat(double x, double y)
    {
        Point pt = new Point(x, y);
        double[] result = _ctFrom.MathTransform.Transform(pt.ToDoubleArray());
        return new Point(result[0], result[1]);
    }

}`

Now Lower48.shp projects to the correct location as it did with methods from SphericalMercator. NY Townships now appear on the map but in the wrong location: near the southern edge of western Africa.
image

Take a look at my custom class and the button code in MainWindow.BtnOSMWithNY_Click().
I am probably overlooking something simple....

pauldendulk commented 5 years ago

The code I posted above already fixed the projection. I now created a pull request with the changes.

It does now show the data at all when zoomed out too much. Not sure what the problem is there.

pauldendulk commented 5 years ago

image

pauldendulk commented 5 years ago

Why do you use render mode wpf?

pauldendulk commented 5 years ago

Why do you use the WPF renderer?

garykindel commented 5 years ago

Using WPF renderer solved a problem I encountered when i tried to create a userControl from the MapSui.ui.wpf.mapcontrol. See Mapsui Error: How to fix PresentationSource is null. An error was being thrown when the view using the usercontrol is created. Purpose behind my usercontrol was to encapsulate any needed code behind the mapsui control. I have since abandoned use of my user control and I am using MapSui.ui.wpf.mapcontrol in my views directly. So short answer:
I don't have a reason anymore why I am using the WPF renderer .

garykindel commented 5 years ago

From your code changes, I can see I had some of my to <=> from between coordinate systems backwards. I also see that my attempt to mimic MapSui.projections.SphericalMercator with Proj.net was less than 100% successful. Thank you for your help. Please feel free to use my code and data as sample code.

pauldendulk commented 5 years ago

ah great. I probably will use your sample. Perhaps I should rewrite the MinimalTransformation to allow adding new transformation entries so users do not need to add a custom transformation.