Accuracy problem

Topics: Algorithms, Data Access, General Topics, SharpMap Project, SharpMap v2.0
Nov 18, 2010 at 4:12 PM

 

I had done an accuracy comparison of converted coordinates (covnerted British National Grid to WGS84) between SharpMap 2 and GeoCoordConversion. I find the GeoCoordConversion produces much better result:

 

Does anyone carry out the similar test? Could anyone share the result?

 

Coordinator
Nov 18, 2010 at 7:51 PM

If you get the SharpMap solution and compile it yourself, you can set configuration to DebugDSProjection or ReleaseDSProjection.

That way you can use another projection library which might produce better results. At least it does for wgs84 -> epsg:3785 (Google Mercator)

Hth FObermaier

Editor
Nov 19, 2010 at 10:25 AM

Sorry to interrupt, but i have some experiance in doing some work on this, but not using SharpMap, so this may not be relant.

What i have found is that most software does not apply a full conversion of the coordinates.

What i mean by that is that the OS suggest a two pass solution. The first pass will perform a basic conversion which places the point in a general vacinity. But that may up to a meter out (depending on where you are located in the UK suprisingly). The second pass takes detailed datum information (which requires a information set from the OS - converted to binary its about 10Mb). The result off applying that shift take the accuracy down to millimeter accuracy (generally).

As an example i wrote a piece of software that would take a google polyline (in WGS94) of a shortest path  and translates it to OSGB36 (BNG). The first pass places the route close to a OS mastermap layer, but not always directly on it (sometimes it was, somethimes it was not). The second pass places it on the road 99% of the time. The remaining 1% is within about 10cm and i suspect that difference is down to differences between the road networks used by google and OS Mastermap.

As i say, having not done any transformation work in Sharpmap (dont even know where to start) so it may not be relavant here.

 

 

 

Coordinator
Nov 19, 2010 at 11:36 AM

If you have code you'd like to share, it is more than welcome

Cheers FObermaier

Editor
Nov 19, 2010 at 4:23 PM

Sure. The system does rely on a large data file (10MB - condensed down from a 43Mb text file that contains same information), thougth.

The solution works with points or limited points. I havent done any large scale tests on it so i'm not sure about performance, allthougth the only difference between a normal transmformation and this should be a X lookup from a table (where X is the number of points in the solution.

It has three classes: One that converts between degrees and radians (it needs all lat and lons to be in radians. But thats basic stuff:

public class Support
    {
        /// <summary>
        /// Converts From Degrees To Radians
        /// </summary>
        /// <param name="degrees">The Value In Degrees</param>
        /// <returns>The Value In Radians</returns>
        public static double DegreeToRadians(double degrees)
        {
            return (degrees * Math.PI) / 180;
        }

        /// <summary>
        /// Converts From Radians To Degrees
        /// </summary>
        /// <param name="radians">The Values In Radians</param>
        /// <returns>The Radians In Degree</returns>
        public static double RadiansToDegrees(double radians)
        {
            return (radians * 180) / Math.PI;
        }
    }
Then theres a second class. Basiclly all that does is hold a look up table containing the shift values from the OS:
    public class Shifts
    {
        public static Single[] shift_e = new Single[876952];
        public static Single[] shift_n = new Single[876952];
        public static Single[] geoid_s = new Single[876952];

        public static void LoadShiftInformation(string fileName)
        {
            string _appPath = fileName;

            BinaryReader reader = new BinaryReader(File.Open(_appPath, FileMode.Open));
            for (int i = 1; i <= 876951; i++)
            {
                Single _east1;
                Single _north1;
                Single _height1;
                _east1 = reader.ReadSingle();
                _north1 = reader.ReadSingle();
                _height1 = reader.ReadSingle();

                shift_e[i] = _east1;
                shift_n[i] = _north1;
                geoid_s[i] = _height1;
            }     
        }
    }
The main trnasformation is like this:
/// <summary>
    /// This class handles the conversion between easting and northings and Lat & Lat (in radians) formats.
    /// </summary>
    public class OSGB36
    {
        public static bool OSGB36Reverse(ref double pLat, ref double pLon, ref double pHeight)
        {
            // Takes a easting, northing and height and converts those locations in the wGs84 (Lat & Lon) as 
            // coorindates

            double eS=0, nS=0;
            double ht_dif=0;
            double E=0, N=0;

            calcShift(pLat, pLon, ref eS, ref nS, ref ht_dif);

            double diff=0, diffn=0, NeweS=0, NewnS=0;
            int iterations = 0;

            do
            {
                E = pLat - eS;         // use shifts to calculate OSGRS80 Grid coordinates
                N = pLon - nS;

                calcShift(E, N, ref NeweS, ref NewnS, ref ht_dif);

                diff = NeweS - eS;
                diffn = NewnS - nS;

                eS = NeweS;
                nS = NewnS;

                iterations = iterations + 1;

            } while ((iterations < 10) && ((Math.Abs(diff) > 0.00001)) || (Math.Abs(diffn) > 0.00001));

            double a = 6378137.0;         // ellipsoid semi minor/major axis
            double b = 6356752.3141;
            double lat0 = 49.0 * (Math.PI / 180.0);   // latitude of true origin (radians)
            double lon0 = -2.0 * (Math.PI / 180.0);   // longitude of true origin (radians)
            double e0 = 400000.0;          // grid Easting of true origin
            double n0 = -100000.0;         // grid Northing of true origin
            double f0 = 0.9996012717;      // scale factor on central meridian

            a *= f0;                         // scale by f0
            b *= f0;

            double n = (a - b) / (a + b);
            double e2 = ((a * a) - (b * b)) / (a * a);
            double latEst = ((N - n0) / a) + lat0;

            double latNew, dLat, sLat, M;
            iterations = 0;

            do
            {
                dLat = latEst - lat0;      // difference in latitude
                sLat = latEst + lat0;      // sum of latitude

                M = ((1 + n + ((5 * n * n) / 4) + ((5 * n * n * n) / 4)) * dLat) -
                       (((3 * n) + (3 * n * n) + ((21 * n * n * n) / 8)) * Math.Sin(dLat) * Math.Cos(sLat)) +
                       ((((15 * n * n) / 8) + ((15 * n * n * n) / 8)) * Math.Sin(2 * dLat) * Math.Cos(2 * sLat)) -
                       (((35 * n * n * n) / 24) * Math.Sin(3 * dLat) * Math.Cos(3 * sLat));
                M *= b;
                latNew = ((N - n0 - M) / a) + latEst;
                latEst = latNew;
                diff = N - n0 - M;
            } while ((iterations++ < 5) && (diff > 0.00001));

            double nu = (a / Math.Pow(1 - e2 * Math.Pow(Math.Sin(latEst), 2.0), 0.5));
            double rho = ((a * (1 - e2)) / Math.Pow(1 - e2 * Math.Pow(Math.Sin(latEst), 2.0), 1.5));
            double n2 = (nu / rho) - 1;

            double VII = Math.Tan(latEst) / (2 * nu * rho);
            double VIII = (Math.Tan(latEst) / (24 * rho * Math.Pow(nu, 3.0))) * (5 + (3 * Math.Pow(Math.Tan(latEst), 2.0)) +
                               n2 - (9 * Math.Pow(Math.Tan(latEst), 2.0) * n2));
            double IX = (Math.Tan(latEst) / (720 * rho * Math.Pow(nu, 5.0))) * (61 + (90 * Math.Pow(Math.Tan(latEst), 2.0)) +
                               (45 * Math.Pow(Math.Tan(latEst), 4.0)));
            double X = 1 / (Math.Cos(latEst) * nu);
            double XI = (1 / (6 * (Math.Cos(latEst) * Math.Pow(nu, 3.0))) * ((nu / rho) + (2 * Math.Pow(Math.Tan(latEst), 2.0))));
            double XII = (1 / (120 * (Math.Cos(latEst) * Math.Pow(nu, 5.0))) * (5 + (28 * Math.Pow(Math.Tan(latEst), 2.0)) + (24 * Math.Pow(Math.Tan(latEst), 4.0))));
            double XIIA = (1 / (5040 * (Math.Cos(latEst) * Math.Pow(nu, 7.0))) * (61 + (662 * Math.Pow(Math.Tan(latEst), 2.0)) +
                               (1320 * Math.Pow(Math.Tan(latEst), 4.0)) + (720 * Math.Pow(Math.Tan(latEst), 6.0))));

            double y = E - e0;
            pLat = latEst - (y * y * VII) + (Math.Pow(y, 4.0) * VIII) - (Math.Pow(y, 6.0) * IX);
            pLon = lon0 + (y * X) - (Math.Pow(y, 3.0) * XI) + (Math.Pow(y, 5.0) * XII) - (Math.Pow(y, 7.0) * XIIA);
            pHeight = pHeight + ht_dif;
            return false;
        }

        public static bool ToOSGB36(ref double pLat,ref double pLon,ref double pHeight)
        {
            double lat=0, lon=0, ht=0;
           
            lat = pLat;
            lon = pLon;
            ht = pHeight;

            double a = 6378137.0;         // ellipsoid semi minor/major axis
            double b = 6356752.3141;
            double lat0 = 49 * (Math.PI / 180.0);   // latitude of true origin (radians)
            double lon0 = -2 * (Math.PI / 180.0);   // longitude of true origin (radians)
            double e0 = 400000.0;          // grid Easting of true origin
            double n0 = -100000.0;         // grid Northing of true origin
            double f0 = 0.9996012717;      // scale factor on central meridian

            a *= f0;                         // scale by f0
            b *= f0;

            double n = (a - b) / (a + b);
            double e2 = ((a * a) - (b * b)) / (a * a);
            double nu = (a / Math.Pow(1 - e2 * Math.Pow(Math.Sin(lat), 2.0), 0.5));
            double rho = ((a * (1 - e2)) / Math.Pow(1 - e2 * Math.Pow(Math.Sin(lat), 2.0), 1.5));
            double n2 = (nu / rho) - 1;

            double dLat = lat - lat0;      // difference in latitude
            double sLat = lat + lat0;      // sum of latitude

            double M;                       // compute Meridian arc (eqn 8.1 [1])
            M = ((1 + n + ((5 * n * n) / 4) + ((5 * n * n * n) / 4)) * dLat) -
                   (((3 * n) + (3 * n * n) + ((21 * n * n * n) / 8)) * Math.Sin(dLat) * Math.Cos(sLat)) +
                   ((((15 * n * n) / 8) + ((15 * n * n * n) / 8)) * Math.Sin(2 * dLat) * Math.Cos(2 * sLat)) -
                   (((35 * n * n * n) / 24) * Math.Sin(3 * dLat) * Math.Cos(3 * sLat));
            M *= b;

            double P = lon - lon0;          // difference in longitude

            // The following are taken from [1] pp10-11
            double I = M + n0;
            double II = (nu / 2) * Math.Sin(lat) * Math.Cos(lat);
            double III = (nu / 24) * Math.Sin(lat) * Math.Pow(Math.Cos(lat), 3.0) * (5 - Math.Pow(Math.Tan(lat), 2.0) + (9 * n2));
            double IIIA = (nu / 720) * Math.Sin(lat) * Math.Pow(Math.Cos(lat), 5.0) *
                               (61 - (58 * Math.Pow(Math.Tan(lat), 2.0)) + Math.Pow(Math.Tan(lat), 4.0));
            double IV = nu * Math.Cos(lat);
            double V = (nu / 6) * Math.Pow(Math.Cos(lat), 3.0) * ((nu / rho) - Math.Pow(Math.Tan(lat), 2.0));
            double VI = (nu / 120) * Math.Pow(Math.Cos(lat), 5.0) * (5 - (18 * Math.Pow(Math.Tan(lat), 2.0)) +
                                Math.Pow(Math.Tan(lat), 4.0) + (14 * n2) - (58 * Math.Pow(Math.Tan(lat), 2.0) * n2));

            double y = I + (Math.Pow(P, 2.0) * II) + (Math.Pow(P, 4.0) * III) + (Math.Pow(P, 6.0) * IIIA);
            double x = e0 + (P * IV) + (Math.Pow(P, 3.0) * V) + (Math.Pow(P, 5.0) * VI);

            double eS,nS,hS;
            eS = 0;
            nS = 0;
            hS = 0;

            calcShift(x, y, ref eS, ref nS,ref hS);

            pLat = x + eS;
            pLon = y + nS;
            pHeight = pHeight - hS;

            return true;
        }

        private static bool calcShift(double x,double y,ref double eS,ref double nS,ref double hS)
        {
            int eIndex = (int)(Math.Floor(x / 1000));
            int nIndex = (int)(Math.Floor(y / 1000));

            // Calculate the values of the shifts at the four corners of the cell

            double se0=0, sn0=0, se1=0, sn1=0, se2=0, sn2=0, se3=0, sn3=0,sg0=0,sg1=0,sg2=0,sg3=0;

            Shift(eIndex, nIndex, ref se0, ref sn0,ref sg0);             // shifts
            Shift((eIndex + 1), nIndex, ref se1, ref sn1,ref sg1);
            Shift((eIndex + 1), (nIndex + 1), ref se2, ref sn2,ref sg2);
            Shift(eIndex, (nIndex + 1), ref se3, ref sn3,ref sg3);

            double X0 = 0, Y0 = 0;

            X0 = eIndex * 1000;
            Y0 = nIndex * 1000;

            double dX = x - X0;
            double dY = y - Y0;
            double t = dX / 1000;
            double u = dY / 1000;

            eS = ((1-t)*(1-u)*(se0)) + ((t)*(1-u)*se1)+((t)*(u)*se2) + ((1-t)*(u)*se3);
            nS = ((1 - t) * (1 - u) * (sn0)) + ((t) * (1 - u) * sn1) + ((t) * (u) * sn2) + ((1 - t) * (u) * sn3);
            hS = ((1 - t) * (1 - u) * (sg0)) + ((t) * (1 - u) * sg1) + ((t) * (u) * sg2) + ((1 - t) * (u) * sg3);
            
            return true;
        }

        /// <summary>
        /// Applies the second level OS shift calculation. Requires the information from
        /// the OSTN02 informaiton file to be loaded.
        /// </summary>
        /// <param name="ei"></param>
        /// <param name="ni"></param>
        /// <param name="eS"></param>
        /// <param name="nS"></param>
        /// <param name="hS"></param>
        /// <returns></returns>
        private static bool Shift(int ei, int ni, ref double eS, ref double nS,ref double hS)
        {
            // Calculate the index in the shift array that will provide the information 
            // required
            int record_number = (ei + (ni * 701) + 1);

             
            // Get the shift infotmation from the look up table.
            eS = (((double)Shifts.shift_e[record_number]));
            nS = ((((double)Shifts.shift_n[record_number])));
            hS = (((double)Shifts.geoid_s[record_number]));
            return true;
        }
     
    }

So to transform a series of points you would do something like this:
OSGB36.Shifts.LoadShiftInformation(path);
OSGB36.OSGB36.ToOSGB36(ref lat,ref lon, ref height);

 

 

 

 

 

 

Editor
Nov 19, 2010 at 4:24 PM

Just to add, theres a detailed description, together with tests on the OS Web Site.