Converting Geographic Coordinates into Lambert2008 coordinates

Question

This is part of a game made in Unity3D.

I wrote this class to convert geographic coordinates into lambert 2008 coordinates based on this reference document:

http://www.ngi.be/Common/Lambert2008/Transformation_Geographic_Lambert_FR.pdf

Reason is that I'm using a Belgian National Geographic Institute map as background, and want to locate waypoint and player position with full accuracy.

Code works fine, and runs in .000001s but I'm unhappy with its readability.

I start by declaring all the constants I need

public const int a = 6378137;
public const double f = 1 / 298.257222101;
public const double lowerLatitude = 49.8333333 * Mathf.Deg2Rad;
public const double upperLatitude = 51.1666667 * Mathf.Deg2Rad;
public const double originLatitude = 50.797815 * Mathf.Deg2Rad;
public const double originLongitude = 4.359215833 * Mathf.Deg2Rad;
public const int originX = 649328;
public const int originY = 665262;

Then I used properties for all values derived from those constants

public static double Excentricity {get => System.Math.Sqrt((2 * f) - (f * f));}
static double MLower { get => M(lowerLatitude);}
static double MUpper { get => M(upperLatitude);}
static double M(double latitude)
{
 return System.Math.Cos(latitude) / System.Math.Sqrt(1 - (Excentricity * Excentricity * System.Math.Pow(System.Math.Sin(latitude), 2)));
}
static double TLower { get => T(lowerLatitude); }
static double TUpper { get => T(upperLatitude); }
static double TOrigin { get => T(originLatitude); }
static double T(double latitude)
{
 return System.Math.Tan(Mathf.PI / 4 - latitude / 2) / System.Math.Pow((1 - Excentricity * System.Math.Sin(latitude)) / (1 + Excentricity * System.Math.Sin(latitude)), Excentricity / 2);
}
static double N { get => (System.Math.Log(MLower) - System.Math.Log(MUpper)) / (System.Math.Log(TLower) - System.Math.Log(TUpper)); }
static double G { get => MLower / (N * System.Math.Pow(TLower, N)); }
static double ROrigin { get => R(TOrigin); }
static double R (double t)
{
 return a * G * System.Math.Pow(t, N);
}

As you can see some of the calculations are written on very long lines. I'm unhappy with their readability, but I have no Idea on how to re-write them to make them look better while not messing up floating point approximations.

Finally I perform the actual conversion using all those calculated values.

public static Vector2 FromGeographicRelative(Vector2 coordinates)
{
 double t = T((double)coordinates.x * Mathf.Deg2Rad);
 double r = R(t);
 double angle = N * ((double)coordinates.y * Mathf.Deg2Rad - originLongitude);
 Vector2 lambertCoord = new Vector2
 {
 x = (float) (r * System.Math.Sin(angle)),
 y = (float) (ROrigin - r * System.Math.Cos(angle))
 };
 return lambertCoord;
}

• jpmc26: There are libraries for coordinate projection. GDAL and ProjNet, for examples. However, .NET is admittedly one of the worst (if not the worst) language for geospatial tools. They tend to be fairly poorly maintained.
• jpmc26: @t3chb0t It's about being practical and not wasting time. A library used more widely than your own code and written by people who deal with this field regularly also has less chance of containing bugs. But if the projection is newer (which it isn't) or if it's buggy in those or doesn't support changing parameters needed (which they do support many), then using it isn't practical.
• user1747281: @jpmc26 Thanks for the info, geospatial localisation is completely new to me, I knew nothing about it 2 days ago, so I didn't know where to look. I very quickly found that reference document with all the formulas I needed to achieve exactly what I wanted, so I didn't even bother looking for a full blown library when I knew I could easily implement it myself. The Belgian National Geographic Institute does provide a windows .dll for conversions, but this has to work on iOS so using their .dll was not an option. I don't like to use a full blown library if I'm only going to use 0.0001% of it
• Eric Duminil: Did you write tests for your conversion? The corresponding EPSG code is 3812, so you can check if your input values are correct here : spatialreference.org/ref/epsg/etrs89-belgian-lambert-2008/html and use this website: epsg.io/transform#s_srs=3812&t_srs=3857 in order to manually create test pairs for your conversion.
• user1747281: I run a test with 5 coordinates: brussels, genk, de panne, bastogne and chimay. I used this tool to determine the target values (uses the official NGI SOAP API). My result is exactly the expected result at every location.

Answer 1

public const double lowerLatitude = 49.8333333 * Mathf.Deg2Rad;

Huh? Mathf.Deg2Rad is a float. Either use Math.PI / 180 or change the type of lowerLatitude to float, but don't work with low-precision values and then implicitly claim that the result is high-precision.

---

static double MLower { get => M(lowerLatitude);}

This (and all the other similar lines) is weird. I can understand using the old-style

static double MLower { get { return M(lowerLatitude); } }

or the new-style

static double MLower => M(lowerLatitude);

but the only reason I can see to use the hybrid form is if there's also a setter.

---

Actually, I think there's a better way to handle derived constants. If you use static readonly fields and a static constructor, the calculation only has to be done once and the methods which you created to share the code can be hidden inside a private scope. Following t3chb0t's comments on variable names:

static readonly double E = System.Math.Sqrt((2 - f) * f);
static readonly double M1;
static readonly double M2;
static readonly double T1;
static readonly double T2;
static readonly double T0;
static readonly double N;
static readonly double G;
static readonly double R0;
static YourClassNameHere()
{
 (double, double) MT(double φ)
 {
 double cosφ = System.Math.Cos(φ);
 double e_sinφ = E * System.Math.Sin(φ);
 double m = cosφ / System.Math.Sqrt(1 - System.Math.Pow(e_sinφ, 2));
 double t = System.Math.Tan(Mathf.PI / 4 - φ / 2) / System.Math.Pow((1 - e_sinφ) / (1 + e_sinφ), E / 2);
 return (m, t);
 }
 double _;
 (_, T0) = MT(φ0);
 (M1, T1) = MT(φ1);
 (M2, T2) = MT(φ2);
 N = System.Math.Log(M1 / M2, T1 / T2);
 G = M1 / (N * System.Math.Pow(T1, N));
 R0 = a * G * System.Math.Pow(T0, N);
}

Note that I've also refactored N to an algebraically equivalent expression (\$\log(x/y)=\log x - \log y\$ and \$\log_b x = \log x / \log b\$) and verified that the result is the same.

---

And as a minor point of spelling: in English the non-circularity of an ellipse is eccentricity with a double-c, reduced from the -xc- of its etymology.

• t3chb0t: I didn't have time to rewrite the formula completely but it looks really nice now, with the original mathematical symbols! One improvement though about (double, double) MT(double φ) - I would give these two names like (double M, double T) etc... or actually it doesn't really matter since you deconstruct them anyway.
• Peter Taylor: @t3chb0t, it's a shame that Unicode subscript numbers aren't also valid in identifiers.
• dfhwze: static readonly seems the best option to me. Simple built-in solution without the need for Lazy, memoization or T4.
• user1747281: You're absolutely right about Deg2Rad, it's legacy from when I initially used UnityEngine.Mathf before switching to System.Math. For the getter syntax, it's the auto-generated code from VS. It generates a property with getter and setter, and I just removed the setter part. Thanks for clarifying this, point taken, lesson learned. Thanks a lot for the refactored code, looks sexy. Looking at it, I can't believe I've not thought of doing it that way from the start. Sometimes the simplest solutions just elude us for some reason.
• user1747281: @PeterTaylor brilliant. Haven't worked much with logs since 12th grade. Do you mind if I edit your answer with that extra bit of explanation?

Answer 2

Since the code is about implementing mathematical formulas I would make an exception and completely ignore the usual naming conventions. Instead, I would use the exact symbols as that document.

This would mean a little bit copy/paste because it's difficult to write symbols like [φ, λ, θ] but the compiler can work with unicode.

I find that code like this:

var rad = 0.01745329252;
var ellipsoide = (a: 6378137, f: 1 / 298.257222101);
var φ1 = 49.8333333 * rad;
var φ2 = 51.1666667 * rad;
var φ0 = 50.797815 * rad;
var λ0 = 4.359215833 * rad;
var x0 = 649328;
var y0 = 665262;

is the most easy one to understand because you can compare it with the book almost without any translation.

• user1747281: Oh, that's awesome, I didn't know the compiler accepted special characters like this. Thanks a lot.
• jpmc26: I would just write phi1 and lambda1. Anyone reasonably familiar with math or physics will know what that means.
• t3chb0t: @jpmc26 hahaha I knew sooner or later you would say that ;-)
• Eric Duminil: I should use unicode more often in my code, it can look awesome. Honest question : Is var a good idea, since OP seems to be mixing floats and doubles?
• t3chb0t: @EricDuminil I think it's ok here becuase double is the default for 0.1 and OP is doing all their calculations with double and only for the final creation of the Vector2 values are casted to float.

Answer 3

• In math heavy code you might use the using static directive, see the microsoft documentation such that System.Math. can be removed throughout the code.
• Excentricity can be made a constant and may have a shorter name as I have seen in other codes dealing with WGS84.
• If possible refer to a document and formula number in a comment
• Although sometime frowned on, I like to layout the the code as much as possible as the mathematical equation looks like. (use extra parentheses if it improves readability)

With this in mind I came to to following example for method T.

using static System.Math;
public const double E2 = (2 * f) - (f * f); # Excentricity squared
public const double E1 = Sqrt(E2); # Excentricity
.....
 static double T(double lat)
 {
 // if possible refer to a document and formula number
 return Tan(PI / 4 - lat / 2) 
 / 
 Pow( ((1 - E1 * Sin(lat)) 
 / 
 (1 + E1 * Sin(lat))) , E1 / 2);
 }

• user1747281: Thanks, a lot for these hints. I'm deliberately not using "using System.Math" because this code is part of a game made with Unity3D. Typically, Unity programmers would use UnityEngine.Mathf, which is a wrapper class for System.Math using floats instead of doubles (because game engines like floats better). So to make it more obvious I explicitly prefix my math calls with System.Math (perhaps a bad practice, would be happy to have your opinion on this) Thanks for the excentricity tip, will definitely do that.
• dfhwze: You can add a using directive: using Math = System.Math
• Jan Kuiken: @user1747281 I am not familiar with Unity programmers conventions, but using static UnityEngine.Mathf might be acceptable, also when using Mathf the constants can be made floats instead of doubles.
• t3chb0t: @user1747281 that this code is part of a game made with Unity3D should have been mentioned in the question so that reviewers have enough context how it is used.

Answer 4

Constants

You have a lot of mathematical constants that aren't constants in your API. Why?

Excentricity, MLower, MUpper, TLower, TUpper, TOrigin, ROrigin, N, G.

EDIT:

In the comments you suggest that the compiler does not accept these constants. But the problem is not the compiler, it's that you should not use the formula to declare a constant. Pre-calculate the constant using T4 and use the result in your derived constants.

// Just paste the pre-calculated result from T4
static double Excentricity = 0.0818191910428158d;

A simpler option is to use static readonlyas Peter suggests in his answer.

Naming Conventions

Some, but not all of your variables and functions are named after mathematical usage. I like it either way, but a combination is an inconsistent naming convention.

public const double f = 1 / 298.257222101;
public const double lowerLatitude = 49.8333333 * Mathf.Deg2Rad;

Aliases

You are using both System.Math and UnityEngine.Mathf. To avoid repeating code I suggest to import the one you used most as a static using and the other as an aliased using.

using static UnityEngine.Mathf;
using Math = System.Math;

• user1747281: Naming conventions: I didn't know you could use mathematical symbols in names, so instead of naming the variable "phi1", "phi2" etc... I opted for a purposeful name. But I now know better thanks to this community :-) Constants: I don't understand what you mean. They're not constants, they're derived from constants. I tried this: public const double Excentricity = System.Math.Sqrt((2 * f) - (f * f)); but the compiler doesn't accept it because Sqrt((2 * f)-(f*f)) is not a constant. Aliases: I like that idea, thanks for the tip
• dfhwze: See my edit => They are constants, you just need to use the formula to get the result and paste that result in your code :) The problem is not the compiler. You should pre-calculate these constants and put the result in the API. For instance, let's take a very easy example: const int x = 1 + 1 -> you could also write const x = 2
• dfhwze: @PeterTaylor I would favor performance optimization over maintainability here. Specially if I feel mathematical constants tend not to change that often (if ever). But I get your point, and I would use a T4 or other template to be able to recalculate these constants whenever I need to. At least use memoization if you keep the formula for a constant value.
• user1747281: Ok I get what you mean, and it's probably the best thing to do in my current situation. However here was my thought process: there's not only one set of lambert conical projection base constants. Here I'm just doing the Lambert2008, but I wanted to keep open the possibility of working with different sets of Lambert constants (Lambert72, french lambert, ....). I could make an excel sheet and calculate excentricity etc whenever I need to do a new set of constants, but it doesn't come across to me as very maintainable code, and if there's a problem with one constant you have to go back to xls
• dfhwze: If you have constants derived from other constants, and you want both performance and maintainability, you should use memoization. Check Lazy<T>.

Answer 5

First thank you all for your help, and specifically @t3chb0t , @JanKuiken , @PeterTaylor and @dfhwze . I learned a lot today.

Here is the refactored code taking into account all the comments. I think it looks a lot better already.

Not only that but performance was multiplied by 20

using static System.Math;
using UnityEngine;
public static class Lambert2008
{
 public const int a = 6378137;
 public const double f = 1 / 298.257222101;
 public const double rad = PI/180;
 public const double φ1 = 49.8333333 * rad;
 public const double φ2 = 51.1666667 * rad;
 public const double φ0 = 50.797815 * rad;
 public const double λ0 = 4.359215833 * rad;
 public const int x0 = 649328;
 public const int y0 = 665262;
 static readonly double E = Sqrt((2 - f) * f);
 static readonly double M1;
 static readonly double M2;
 static readonly double T1;
 static readonly double T2;
 static readonly double T0;
 static readonly double N;
 static readonly double G;
 static readonly double R0;
 static Lambert2008()
 {
 double _;
 (_, T0) = MT(φ0);
 (M1, T1) = MT(φ1);
 (M2, T2) = MT(φ2);
 N = Log(M1 / M2, T1 / T2);
 G = M1 / (N * Pow(T1, N));
 R0 = a * G * Pow(T0, N);
 }
 static (double, double) MT(double φ)
 {
 var e_sinφ = E * Sin(φ);
 var m = Cos(φ) / Sqrt(1 - Pow(e_sinφ, 2));
 var t = CalculateT(φ);
 return (m, t);
 }
 static double R(double t)
 {
 return a * G * Pow(t, N);
 }
 static double CalculateT(double φ)
 {
 var e_sinφ = E * Sin(φ);
 return Tan(PI / 4 - φ / 2)
 / Pow((1 - e_sinφ) / (1 + e_sinφ), E / 2);
 }
 public static Vector2 FromGeographicRelative(Vector2 coordinates)
 {
 double _, t, r;
 (_, t) = MT(coordinates.x * rad);
 r = R(t);
 var θ = N * (coordinates.y * rad - λ0);
 var res = new Vector2
 {
 x = (float)(r * Sin(θ)),
 y = (float)(R0 - r * Cos(θ))
 };
 return res;
 }
 public static Vector2 FromGeographic(Vector2 coordinates)
 {
 return FromGeographicRelative(coordinates) + new Vector2(x0, y0);
 }
}

• t3chb0t: Coool! This has become really nice!
• user1747281: Also, I just tested performance. 10000 calls to FromGeographicRelative as benchmark. Before refactoring: 0.09s; After refactoring: 0.005s so about 20x performance gain
• Peter Taylor: Might as well refactor MT back to just M. The only reason I combined them in one was to save recalculating e_sinφ.