Преобразовать координаты OSGB 36 в широту / долготу

Question

Я хочу преобразовать британские OSGB 36 координаты в WGS 84 (т.е. "стандартные" широта и долгота), чтобы отобразить их в файл KML для Google Earth .

Как лучше всего это сделать? Я внедряю в VB .NET.

Я, вероятно, должен добавить, что мой вопрос не «Как мне написать файл KML?». Мой вопрос: «Как я могу конвертировать между этими двумя системами координат?» !!

Я надеялся, что найдется библиотека, которую я мог бы использовать, вместо того, чтобы использовать свою собственную функцию - похоже, что-то вроде того, что кто-то другой мог бы реализовать.

Answer 1

Компания, в которой я работаю, только что открыла библиотеку, чтобы сделать именно это: http://code.google.com/p/geocoordconversion/

Answer 2

//=======================================================================
/* Project: Geocode.Service
*  Author: Steve Loughran
*  Copyright (C) 2000 Steve Loughran
*  See license.txt or license.html for license and copyright information 
*  RCS $Header: /cvsroot/geocode/geocode.net/src/library/Osgb.cs,v 1.4 2002/04/23 05:12:37 steve_l Exp $
*  jedit:mode=csharp:tabSize=4:indentSize=4:syntax=true:
*/
//=======================================================================


using System;

namespace net.sourceforge.geocode {

/// <summary>
/// OSGB conversion routines. xlated from C++ to Java to C#
/// </summary>

public class OsgbGridReference: GeoMath
{

private string _gridSquare="";
private long _easting=0;
private long _northing=0;

/// <summary>
/// empty constructor
/// </summary>
public OsgbGridReference() {}

/// <summary>
/// constructor from gridref
/// </summary>
public OsgbGridReference(String gridSquare, 
  long easting,
  long northing) {
 SetGridRef(gridSquare,northing,easting);
}

/// <summary>
/// constructor from co-ordinates
/// </summary>
public OsgbGridReference(double latitude, double longitude) {
 SetPosition(latitude,longitude);
}

/// <summary>
/// constructor from position
/// </summary>
public OsgbGridReference(Position position)
 : this(position.Latitude,position.Longitude) {
}

/// <summary>grid square property</summary>    
public string GridSquare {
 get {return _gridSquare;}
 set {_gridSquare=value;}
}
/// <summary>northing property</summary>    

public long Northing {
 get {return _northing;}
 set {_northing=value;} 
}

/// <summary>easting property</summary>    
public long Easting {
 get {return _easting;}
 set {_easting=value;} 
}

/// <summary>
/// set the grid refernce
/// </summary>
/// <returns> </returns>

public void SetGridRef(String gridSquare, 
  long northing, 
  long easting) {
 _gridSquare=gridSquare;
 _northing=northing;
 _easting=easting;
}

/// <summary>
///  rudimentary validity test. A better one is to
/// extract the position
/// </summary>
/// <returns>true iff there is a gridsquare </returns>

public bool IsValid() 
 {return _gridSquare.Length==2;}


/// <summary>
/// get a position object from a location
/// </summary>
/// <returns> Position </returns>

public Position ToPosition() {
 double lat,lon;
 ConvertOSGBtoLL(_gridSquare,_northing,_easting,
   out lat, out lon);
 Position p=new Position(lat,lon);
 p.Name=ToString();
 return p;
}   

/// <summary>
/// set a gridref from a lat/long tuple
/// </summary>
/// <returns> success flag </returns>

public bool SetPosition(double latitude, double longitude) {
 _gridSquare=ConvertLLtoOSGB(latitude,
   longitude,
   out _northing,
   out _easting);
 return IsValid();
}

/// <summary>
/// set a gridref from a position
/// </summary>
/// <returns> success flag </returns>

public bool SetPosition(Position position) {
 return SetPosition(position.Latitude,position.Longitude);
}

/// <summary>
///  stringify
/// </summary>

public override string ToString() {
 return String.Format("{0} {1:000} {2:000}",
  _gridSquare,Easting,Northing);
}

/// <summary>
///  equality test: works on lat and long
/// </summary>

public override bool Equals(object o) {
 OsgbGridReference pos=(OsgbGridReference)o;
 return _gridSquare==pos._gridSquare &&
  _northing==pos._northing &&
  _easting==pos._easting;
}

/// <summary>
/// hash code builder
/// </summary>

public override int GetHashCode() {
        return (int)(_easting+_northing); 
}

/// <summary>
/// determines base co-ordinates of a square like "ST"
/// </summary>
///    <parameter name="OSGBGridSquare"> square </parameter>
///    <parameter name="easting"> easting</parameter>  
///    <parameter name="northing"> northing</parameter> 
/// <returns>true if the coordinates were in range</returns>

static bool ConvertOSGBSquareToRefCoords(string OSGBGridSquare,
         out long  easting,
                           out long  northing) {
 int pos, x_multiplier=0, y_multiplier=0;
 string GridSquare = "VWXYZQRSTULMNOPFGHJKABCDE";
 bool trouble=false;
 long east,north;
 easting=northing=0;
 //find 500km offset
    OSGBGridSquare=OSGBGridSquare.ToUpper();
 char ch = OSGBGridSquare[0];
 switch(ch) {
  case 'S': x_multiplier = 0; y_multiplier = 0; break;
  case 'T': x_multiplier = 1; y_multiplier = 0; break;
  case 'N': x_multiplier = 0; y_multiplier = 1; break;
  case 'O': x_multiplier = 1; y_multiplier = 1; break;
  case 'H': x_multiplier = 0; y_multiplier = 2; break;
  case 'J': x_multiplier = 1; y_multiplier = 2; break;
        default: trouble=true; break;
 }
    if(!trouble) {
  east=x_multiplier * 500000L;
  north=y_multiplier * 500000L;
  //find 100km offset and add to 500km offset to get coordinate of
  //square point is in
  char subsquare=OSGBGridSquare[1];
  pos = GridSquare.IndexOf(subsquare);
     if(pos>-1) {
   east += ((pos % 5) * 100000L);
   north += ((pos / 5) * 100000L);
     easting=east;
   northing=north;
     }
     else {
         trouble=true;
     }
    }
    return !trouble;
}

///<summary>
///convert a internal OSGB coord to lat/long
///Equations from USGS Bulletin 1532
///East Longitudes are positive, West longitudes are negative.
///North latitudes are positive, South latitudes are negative
///Lat and Long are in decimal degrees.
///Written by Chuck Gantz- chuck.gantz@globalstar.com
///</summary>
///    <parameter name="OSGBEasting">easting </parameter> 
///   <parameter name="OSGBEasting">northing </parameter> 
///   <parameter name="OSGBZone"> gridsquare</parameter>  
///   <parameter name="latitude"> latitude</parameter> 
///   <parameter name="longitude"> longitude</parameter> 

static void ConvertOSGBtoLL(string OSGBZone,
   double OSGBNorthing,
   double OSGBEasting,
   out double latitude,
   out double longitude) {
 double k0 = 0.9996012717;
 double a;
 double eccPrimeSquared;
 double N1, T1, C1, R1, D, M;
 double LongOrigin = -2;
 double LatOrigin = 49;
 double LatOriginRad = LatOrigin * deg2rad;
 double mu, phi1, phi1Rad;
 double x, y;
    long northing;
    long easting;

 //Airy model of the ellipsoid.
 double majoraxis = a = 6377563.396;
 double minoraxis = 6356256.91;
 double eccSquared = (majoraxis * majoraxis - minoraxis * minoraxis) /
       (majoraxis * majoraxis);
 double e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
 //only calculate M0 once since it is based on the origin of the OSGB projection, which is fixed
 double  M0 = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatOriginRad
    - (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatOriginRad)
    + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatOriginRad)
    - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatOriginRad));
 ConvertOSGBSquareToRefCoords(OSGBZone, out easting, out northing);
 //remove 400,000 meter false easting for longitude
 x = OSGBEasting - 400000.0 + easting;
 //remove 100,000 meter false easting for longitude
 y = OSGBNorthing + 100000.0 + northing;
 eccPrimeSquared = (eccSquared)/(1-eccSquared);
 M = M0 + y / k0;
 mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/
      64-5*eccSquared*eccSquared*eccSquared/256));
 phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
    + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
    +(151*e1*e1*e1/96)*sin(6*mu);
 phi1 = phi1Rad*rad2deg;
 N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad));
 T1 = tan(phi1Rad)*tan(phi1Rad);
 C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad);
 R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5);
 D = x/(N1*k0);
 latitude = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
     +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
 latitude *= rad2deg;
 longitude = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
     *D*D*D*D*D/120)/cos(phi1Rad);
 longitude = LongOrigin + longitude * rad2deg;
}


/// <summary>
/// converts lat/long to OSGB coords.  Equations from USGS Bulletin 1532
/// East Longitudes are positive, West longitudes are negative.
/// North latitudes are positive, South latitudes are negative
/// Lat and Long are in decimal degrees
/// </summary>
///  Written by Chuck Gantz- chuck.gantz@globalstar.com
///   <parameter name="latitude"> IN latitude</parameter> 
///   <parameter name="longitude">IN longitude </parameter> 
///   <parameter name="OSGBEasting"> OUT easting</parameter> 
///  <parameter name="OSGBNorthing"> OUT northing</parameter> 

static public string ConvertLLtoOSGB(double latitude,
    double longitude,                
                out long OSGBNorthing,
    out long OSGBEasting) {
 double a;
 double eccSquared;
 double k0 = 0.9996012717;
 double LongOrigin = -2;
 double LongOriginRad = LongOrigin * deg2rad;
 double LatOrigin = 49;
 double LatOriginRad = LatOrigin * deg2rad;
 double eccPrimeSquared;
 double N, T, C, A, M;
 double LatRad = latitude*deg2rad;
 double LongRad = longitude*deg2rad;
 double easting, northing;
 double majoraxis = a = 6377563.396;//Airy
 double minoraxis = 6356256.91;//Airy
 eccSquared = (majoraxis * majoraxis - minoraxis * minoraxis) /
       (majoraxis * majoraxis);
 //only calculate M0 once since it is based on the origin
 //of the OSGB projection, which is fixed
 double  M0 = a*((1  - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatOriginRad
    - (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatOriginRad)
    + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatOriginRad)
    - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatOriginRad));
 eccPrimeSquared = (eccSquared)/(1-eccSquared);
 N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad));
 T = tan(LatRad)*tan(LatRad);
 C = eccPrimeSquared*cos(LatRad)*cos(LatRad);
 A = cos(LatRad)*(LongRad-LongOriginRad);
 M = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64- 5*eccSquared*eccSquared*eccSquared/256)*LatRad
    - (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
    + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
    - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad));
 easting = (double)(k0*N*(A+(1-T+C)*A*A*A/6
    + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120));
 easting += 400000.0; //false easting
 northing = (double)(k0*(M-M0+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
     + (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720)));
 northing -= 100000.0;//false northing
 return ConvertCoordsToOSGBSquare(easting, northing,out OSGBEasting, out OSGBNorthing);
}


/// <summary>
/// convert a internal OSGB coord to a gridsquare and internal values.
/// </summary>
///    <parameter name="easting"> easting</parameter> 
///    <parameter name="northing"> northing</parameter>
///    <parameter name="OSGBEasting">OSGBEasting out </parameter>
///    <parameter name="OSGBNorthing">OSGBNorthing out </parameter>

static string ConvertCoordsToOSGBSquare(double easting,
     double northing,
                    out long  OSGBEasting,
                    out long  OSGBNorthing)
{
 string GridSquare = "VWXYZQRSTULMNOPFGHJKABCDE";
 long posx, posy; //positions in grid
 OSGBEasting = (long)(easting + 0.5); //round to nearest int
 OSGBNorthing = (long)(northing + 0.5); //round to nearest int
 string OSGBGridSquare="";

 //find correct 500km square
 posx = OSGBEasting / 500000L;
 posy = OSGBNorthing / 500000L;
 if(posx<0 || posx>4 || posy<0 || posy>4)
  return "";
 long offset=posx + posy * 5 + 7;
    OSGBGridSquare+= GridSquare[(int)offset];

 //find correct 100km square
 posx = OSGBEasting % 500000L;//remove 500 km square
 posy = OSGBNorthing % 500000L;//remove 500 km square
 posx = posx / 100000L;//find 100 km square
 posy = posy / 100000L;//find 100 km square
 if(posx<0 || posx>4 || posy<0 || posy>4)
  return "";
 offset=posx + posy * 5; 
 OSGBGridSquare+= GridSquare[(int)offset];

 //remainder is northing and easting
 OSGBNorthing = OSGBNorthing % 500000L;
 OSGBNorthing = OSGBNorthing % 100000L;
 OSGBEasting = OSGBEasting % 500000L;
 OSGBEasting = OSGBEasting % 100000L;
    return OSGBGridSquare;
}



/// <summary>
/// turn the latitude and longitude into a string
/// </summary>
///    <parameter name="latitude"> lat</parameter>
///    <parameter name="longitude"> long</parameter> 
///    <parameter name="infill"> text between coordinates</parameter>  
///    <returns>return something like E004 N123</returns>

static string GetSensibleLatLongstring(double latitude,
  double longitude,
        int decimals,
        string infill) {
    bool bNorth=latitude>=0;
 bool bWest=longitude<=0;
    latitude=Math.Abs(latitude);
    longitude=Math.Abs(longitude);
    double multiplier=(int)pow(10,decimals);
 int hiLat=(int)latitude;
    int lowLat=(int)((latitude-hiLat)*multiplier);
    double newLat=lowLat/multiplier+hiLat;
 int hiLong=(int)longitude;
    int lowLong=(int)((longitude-hiLong)*multiplier);
    double newLong=lowLong/multiplier+hiLong;
 return (bNorth?"N":"S")+newLat+infill+
  (bWest?"W":"E")+newLong;
}




/* legacy java test harness
  public static void main(string[] args)
 {
    string message;
    if(args.length!=3)
     {
        message="need a grid reference like ST 767 870";
        }
    else
     {
        LongRef north=new LongRef();
        LongRef east=new LongRef();
        bool b=ConvertOSGBSquareToRefCoords(args[0],east,north);
        double easting=Double.valueOf(args[1]).doubleValue();
        double northing=Double.valueOf(args[2]).doubleValue();
        DoubleRef rlatitude=new DoubleRef ();
        DoubleRef rlongitude=new DoubleRef ();
        OSGBtoLL(easting,northing,args[0],rlatitude,rlongitude);
        double latitude=rlatitude.getValue();
        double longitude=rlongitude.getValue();

        bool bNorth=latitude>=0;
        bool bWest=longitude<=0;

  message="Latitude & Longitude ="+latitude+" , " + longitude;
        message+="\r\n or "+GetSensibleLatLongstring(latitude,
             longitude,
                            3,
                            " ");
  string square=new string();
  square=LLtoOSGB(latitude,longitude,north,east);
  message+="\r\n Grid ="+square+" "+east+" "+north;
//       message="evaluation failure on "+args[0];
        }
        System.out.print(message);
    }
*/


}; //class
};//namespace

Answer 3

Во-первых, согласно этой странице, связанной с OSGB 36 :

Миф 4: «Есть точные математические формулы для изменения между системами координат»

Во-вторых, по той же ссылке: « Из одной системы координат в другую: геодезические преобразования » включает в себя раздел «Приблизительное преобразование WGS84 в OSGB36 / ODN»

РЕДАКТИРОВАТЬ: Обратите внимание, когда ОС говорит "приблизительный", они означают с ошибками> 5 м.

Answer 4

Вам необходимо реализовать преобразование Гельмерта . Я написал конвертацию в Javascript, которую вы можете адаптировать. Алгоритм, используемый сценарием для преобразований WGS84-OSGB36, получен из таблицы OSGB с разрешения. Точность преобразования составляет порядка 7 м для 90% территории Великобритании и должна быть такой же, как и при обычном приемнике GPS.

См. документацию и источник для получения более подробной информации.

Изменить: вы можете проверить это OCX , который включает в себя источник.

Answer 5

Мы используем библиотеку proj.4 для преобразования координат WGS84 <-> OSGB36 <-> OSGBGRID, и она работает очень хорошо. Но мы используем C ++, поэтому я не знаю, сможете ли вы заставить его работать под VB.NET. Там могут быть обертки или что-то? (На ссылке выше упоминается PROJ.4 VB Wrappers).