//=======================================================================
/* 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
