Add offset to WGS84 coords

Topics: Algorithms, SharpMap Project
Sep 11, 2012 at 2:04 PM

Hi all, I have this problem:

I have some WGS84 coords and I have to add an offset to these coords, for example:

lat: 

12,3913578190804

lng:

43,9950902282485

Is there any method to add for example 10 meters to lat and 5 meters to lng?

By now I use this code:

1) Using dotspatial I convert the coords to Gauss-Boaga (some random decimal coordinate system)

2) Add the offset to the converted coords

3) Using dotspatial I reconvert the coords from Gauss-Boaga to WGS84

With this, I get a 3mm/m error

Coordinator
Sep 11, 2012 at 2:42 PM

How about a pure mathematica approach:

// Reference:
// http://projnet.codeplex.com/discussions/77458
// Post by irezax
double[] ToMercator(double lon, double lat)
{
    double x = lon * 20037508.34 / 180;
    double y = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180);
    y = y * 20037508.34 / 180;
    return new double[] {x, y};
}

double[] ToLonLat(double x, double y) 
{
    double lon = (x / 20037508.34) * 180;
    double lat = (y / 20037508.34) * 180;

    lat = 180/Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2);
    return new double[] {lon, lat};
}

Hth FObermaier

Sep 12, 2012 at 8:09 AM

So, adding offset to my coords would be like this (assuming +10m lon, +0m lat)

 

double[] xy = new double[] { 12.3798088936852, 43.993187787765493 };


double lon = (10 / 20037508.34) * 180;                        
double lat = (0 / 20037508.34) * 180;
lat = 180 / Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2);                        
double[] xy1 = new double[] { xy[0] + lon, xy[1] + lat };

so my point become 12.379898725213623,43.993187787765493
I calculated the distance from my point and this point with Haversine formula, and the result is 7.1926152788263975 meters

 

If I do this (+0m lon, +10m lat)
double lon = (0 / 20037508.34) * 180;                        
double lat = (10 / 20037508.34) * 180;
lat = 180 / Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2);                        
double[] xy1 = new double[] { xy[0] + lon, xy[1] + lat };
my point become 12.3798088936852,43.993277619293913
and the distance is 9.9986380096607643 (1,3mm per meter error)


Sep 12, 2012 at 8:37 AM

I tried this also in google maps  ... I tried to move my geometries 10m to north and then 10m to west, and the latitude offset is quite perfect, but the longitude offset is ~7m

Coordinator
Sep 12, 2012 at 9:25 AM
Edited Sep 12, 2012 at 10:34 AM

Does the follow up by trieuvy in the above mentioned post help:

 

double[] ToMercator2(double lon, double lat)
{
    double x = lon * 20037508.34 / 180;
    double y = -180d;
    if (lat > -90)
    y = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180);
    y = y * 20037508.34 / 180; return new double[] { x, y };
}

Edit:

If not you may find this useful:

  • http://williams.best.vwh.net/avform.htm#LL
  • http://williams.best.vwh.net/gccalc.htm

Hth FObermaier

Developer
Sep 14, 2012 at 3:41 AM

 Hi steelraiden,

20037508.34 = Math.PI * (radius of the earth in meter).
You need calculate radius of the earth in mm axactly.
B/c: replace
20037508.34 = Math.PI*6378100 ( reasonable adjust the radius of the earth to get accurate results that you need).

TrieuVy.
Sep 18, 2012 at 12:23 PM

Hi trieuvy,

I didn't understand where I have to write your code... Maybe here?


from

double lon = (0 / 20037508.34) * 180;                        
double lat = (10 / 20037508.34) * 180;
lat = 180 / Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2);                        
double[] xy1 = new double[] { xy[0] + lon, xy[1] + lat };


to


double lon = (0 / Math.PI * 6378100) * 180;                        
double lat = (10 / Math.PI * 6378100) * 180;
lat = 180 / Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2);                        
double[] xy1 = new double[] { xy[0] + lon, xy[1] + lat };

Developer
Sep 20, 2012 at 1:31 AM
Edited Sep 20, 2012 at 1:34 AM

Yes. try some thing like that.

this functions using spherical trigonometry and assume on a perfect sphere.

// Reference:
// http://projnet.codeplex.com/discussions/77458
// Post by irezax
double[] ToMercator(double lon, double lat)
{
    double x = lon * 20037508.34 / 180;
    double y = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180);
    y = y * 20037508.34 / 180;
    return new double[] {x, y};
}

double[] ToLonLat(double x, double y) 
{
    double lon = (x / 20037508.34) * 180;
    double lat = (y / 20037508.34) * 180;

    lat = 180/Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2);
    return new double[] {lon, lat};
}

Reference: http://alastaira.wordpress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection/

Edited: do you want to calculate the latitude no on a perfect sphere?

Oct 8, 2012 at 3:27 PM

I try to better explain my problem...

I have to convert a layer to .KML.
Only WGS84 coords are allowed in .KML file.
So I use DotSpatial to reproject my coords to WGS84

Now: sometimes the conversion is not good because of the source, and my polygons appear far from the point that they have to be.

To correct the error, I thought to apply an offset in the conversion:

//for all vertices in all polygons
double[] xy = new double[] { myX+offsetx, myY+offsetY };
Reproject.ReprojectPoints(xy, new double[0], _OldCoords, _NewCoords, 0, 1);


now: the offset has to be in the same coordinate system of my point, so it has to be in WGS84. 
My program is used by some people that don't understand the difference from a coordinate system to another, so I thought to ask them an offset in meters.
After that, I thought convert my offset in WGS84 and add it to my coords


Is that possible?

Coordinator
Oct 9, 2012 at 7:46 AM
steelraiden wrote:

 Now: sometimes the conversion is not good because of the source, and my polygons appear far from the point that they have to be.  

If you have layers that do not have the conversion problem when creating the WGS84 KML coordinates, maybe assigning another projection to the layers that don't match would be a better idea. On the Proj.Net site, there is an elaborate post on how to setup the wkt for EPSG:3857 to make it work correctly.