Search code examples
c#map-projections

Take latitude/longitude on a Lambert Conic Conformal Projection Map and get pixel X/Y

I am implementing a map in .net 4.5 WPF application.

I was given an image of a map in lambert conic conformal projection.

I am struggling with the conversion from lat/long to lambert conic so that I can plot points on my map.

My image is a size of 3960x2452 and will be used as the background of a Canvas. I am trying to take lat/long from datapoints and plot them on this canvas. Which is why I need to convert to lambert conic/pixel xy

public LambertConicConversions(double originLatitude, double originLongitude, double standardParallel1, double standardParallel2, double centralMeridian) 
{ 
    this.originLatitude = this.DegreesToRadians(originLatitude);
    this.originLongitude = this.DegreesToRadians(originLongitude);
    this.standardParallel1 = this.DegreesToRadians(standardParallel1);
    this.standardParallel2 = this.DegreesToRadians(standardParallel2);
    this.centralMeridian = this.DegreesToRadians(centralMeridian);

    double m1 = Math.Cos(this.standardParallel1) / Math.Sqrt(1 - eccentricity * Math.Pow(Math.Sin(this.standardParallel1), 2));

    double m2 = Math.Cos(this.standardParallel2) / Math.Sqrt(1 - eccentricity * Math.Pow(Math.Sin(this.standardParallel2), 2));

    double t1 = Math.Tan(Math.PI / 4 - this.standardParallel1 / 2) / 
        Math.Pow(1 - Math.Sqrt(eccentricity) * Math.Sin(this.standardParallel1) / 
        (1 + Math.Sqrt(eccentricity) * Math.Sin(this.standardParallel1)), Math.Sqrt(eccentricity) / 2);

    double t2 = Math.Tan(Math.PI / 4 - this.standardParallel2 / 2) / 
        Math.Pow((1 - Math.Sqrt(0.0818191908426) * Math.Sin(this.standardParallel2)) / 
        (1 + Math.Sqrt(0.0818191908426) * Math.Sin(this.standardParallel2)), Math.Sqrt(0.0818191908426) / 2);

    double t0 = Math.Tan(Math.PI / 4 - DegreesToRadians(this.originLatitude) / 2) / 
        Math.Pow((1 - Math.Sqrt(0.0818191908426) * Math.Sin(DegreesToRadians(this.originLatitude))) / 
        (1 + Math.Sqrt(0.0818191908426) * Math.Sin(DegreesToRadians(this.originLatitude))), Math.Sqrt(0.0818191908426) / 2);

    n = (Math.Log(m1) - Math.Log(m2)) / (Math.Log(t1)- Math.Log(t2));
    F = m1 / (n * Math.Pow(t1, n));
    this.rho0 = earthRadius * F * Math.Pow(t0, n);
}

public void LatLonToLambertXY(double latitude, double longitude, out double x, out double y)
{
    latitude = this.DegreesToRadians(latitude);
    longitude = this.DegreesToRadians(longitude);

    double t = 1 / Math.Tan(Math.PI / 4 - latitude / 2) / Math.Pow((1 - Math.Sqrt(eccentricity) * Math.Sin(latitude)) / (1 + Math.Sqrt(eccentricity) * Math.Sin(latitude)), Math.Sqrt(eccentricity) / 2);

    rho = earthRadius * F * Math.Pow(t, n);
    double theta = n * (longitude - originLongitude);

    double lambertX = rho * Math.Sin(theta);
    double lambertY = rho0 - rho * Math.Cos(theta);

    x = lambertX;
    y = lambertY;
}

Here is the setup code for my conversions. Standard parallels of 33 and 45 Central Meridian of -90 Origin Latitude of 39 for center. Map spans from -130 to -65 long and 50 to 25 latitude. Preferably my origin be 0,0

I have tried using ProjNet3GeoApi and it gives really weird results. I have used the formulas available regarding lambert conic.

When I run these my x coordinates seem relatively close and increment correctly.

My y values start opposite and go from bottom up instead of 0 down.

Solution

  • Not sure about the code in your question, but general map projections (e.g. those provided by the ProjNet4GeoAPI library) transform from geodetic coordinates (latitude/longitude) to cartesian coordinates and back.

    The units of the cartesian coordinate system are typically meters measured from the origin of the projection, where positive x value go to the right and positive y values go upwards.

    A map projection does not transform from meters to pixels of a computer display. This is something that has to be done afterwards and needs a view scale (from meters to pixels) and a view origin (in meters or in pixels) as input. By using a negative value for scaling the y-axis, you would invert the direction of the y coordinates.

    With using ProjNet4GeoAPI, a transform class could look like shown below. Besides the IMathTransform that transforms geodetic to cartesian coordinates it has three properties for the view transform and a method LocationToView that combines both. The view origin is given in meters (from the upper left corner of the map viewport to the origin of the projected coordinate system).

    Also note that the values of the projection parameters are all hardcoded in the WKT string. They would certainly be provided as constructor parameters in a real-world application. You must however be aware that setting False_Easting and False_Northing only shifts the origin of the projected coordinate system. It does not invert the direction of its axes, i.e. you can not make y values go downwards by setting a False_Northing.

    public class LambertConicalTransform
{
    public LambertConicalTransform()
    {
        var wkt
            = "PROJCS[\"Lambert_Conformal_Conic\","
            + "GEOGCS[\"GCS_WGS_1984\","
            + "DATUM[\"WGS_1984\","
            + "SPHEROID[\"WGS_84\",6378137.0,298.257223563],"
            + "TOWGS84[0,0,0,0,0,0,0]],"
            + "PRIMEM[\"Greenwich\",0.0],"
            + "UNIT[\"Degree\",0.0174532925199433]],"
            + "PROJECTION[\"Lambert_Conformal_Conic_2SP\"],"
            + "PARAMETER[\"False_Easting\",0.0],"
            + "PARAMETER[\"False_Northing\",0.0],"
            + "PARAMETER[\"Central_Meridian\",-90],"
            + "PARAMETER[\"Standard_Parallel_1\",33],"
            + "PARAMETER[\"Standard_Parallel_2\",45],"
            + "PARAMETER[\"Latitude_Of_Origin\",39],"
            + "UNIT[\"Meter\",1.0]]";

        var cs = new CoordinateSystemFactory().CreateFromWkt(wkt);

        LocationToMapTransform = new CoordinateTransformationFactory()
            .CreateFromCoordinateSystems(GeographicCoordinateSystem.WGS84, cs)
            .MathTransform;
    }

    public IMathTransform LocationToMapTransform { get; }
    public double ViewScale { get; set; } = 1;
    public double ViewOriginX { get; set; }
    public double ViewOriginY { get; set; }

    public Coordinate LocationToView(double longitude, double latitude)
    {
        return LocationToView(new Coordinate(longitude, latitude));
    }

    public Coordinate LocationToView(Coordinate location)
    {
        var mapPoint = LocationToMapTransform.Transform(location);

        return new Coordinate(
            ViewScale * (ViewOriginX + mapPoint.X),
            ViewScale * (ViewOriginY - mapPoint.Y)); // inverts y direction
    }
}

    Now let's use this transform class to produce some sample data with a view scale of 1e-4 (px/m) and a view origin at (x=4500 km, y=2000 km), in a simple UI application (WPF here).

    In your application, you have to find out about the scale of your map image in order to use the correct value for the view scale. The image's width and height and alignment in your UI determines the values of the origin x and y.

    The view model:

    public class ViewModel
{
    public List<Point> Points { get; } = new();

    public LambertConicalTransform Transform { get; } = new()
    {
        ViewScale = 1e-4,
        ViewOriginX = 4500000,
        ViewOriginY = 2000000,
    };

    public void CreatePoints()
    {
        for (double lon = -130; lon <= -65; lon += 1)
        {
            var p1 = Transform.LocationToView(new Coordinate(lon, 50));
            var p2 = Transform.LocationToView(new Coordinate(lon, 25));
            Points.Add(new Point(p1.X, p1.Y));
            Points.Add(new Point(p2.X, p2.Y));
        }

        for (double lat = 25; lat <= 50; lat += 1)
        {
            var p1 = Transform.LocationToView(new Coordinate(-130, lat));
            var p2 = Transform.LocationToView(new Coordinate(-65, lat));
            Points.Add(new Point(p1.X, p1.Y));
            Points.Add(new Point(p2.X, p2.Y));
        }
    }
}

    The view's code behind:

    public partial class MainWindow : Window
{
    public MainWindow()
    {
        InitializeComponent();

        var vm = new ViewModel();
        vm.CreatePoints();
        DataContext = vm;
    }
}

    And finally the view:

    <ItemsControl ItemsSource="{Binding Points}">
    <ItemsControl.ItemsPanel>
        <ItemsPanelTemplate>
            <Canvas/>
        </ItemsPanelTemplate>
    </ItemsControl.ItemsPanel>
    <ItemsControl.ItemTemplate>
        <DataTemplate>
            <Path Fill="Black">
                <Path.Data>
                    <EllipseGeometry
                        Center="{Binding}"
                        RadiusX="1.5" RadiusY="1.5"/>
                </Path.Data>
            </Path>
        </DataTemplate>
    </ItemsControl.ItemTemplate>
</ItemsControl>

    The produced result:

    enter image description here