The Maidenhead System

The maidenhead locator was invented to provide a worldwide simple to use locator based on latitude and longitude.

The format is recursive with alternating pairs of numbers and letters. The first character representing Latitude and the second Longitude. Resolution increases with the number of pairs of digits. More than 4 pairs is unusual.

The picture below generated by the excellent AZ_Proj utility gives the idea.

The OSGB system

The UK NGR is based on a projection of an ellipsoid to allow a curved surface to be represented on a flat map. It is only valid for the UK.

The ellipsoid used in the OSGB36 datum which was invented in 1830 by the Astronomer Royal, Sir George Airy. The caluclations in this programme are fairly detailed and should give results accurate to a few metres. Beware of Ellipsoids. GPS data is often based on WGS84. Using the incorrect ellipsoid can result in errors of several 100m.

The UK national grid reference should be written as a two letter plus 4 digit, 6 digit or 8 digit locator.

E.g. SU59, SU500900, SU50009000

The more digits the higher the accuracy - 8 is generally more than enough resolution unless you are into surveying. These locators are commonly printed as a grid on UK Ordanance Survey maps.

The grid covering the UK is 700km by 1300km. The lettters indicate the large square, the first letter represents a 500x500km square and can be S,T, N, O, H or J. Placed over the UK as below, with North being at the top.

H J
N O
S T

The second letter represents a 100kmx100km square within the 500km squares. It can be any letter except I. The format is:

A B C D E
F G H J K
L M N O P
Q R S T U
V W X Y Z

The numbers are then the range East and North of the bottom left corner of the square. With 2 digits, the first is East and the next is North. Naturally, the digit is in units of 10km. With 2 digits, the first 2 are East and the next 2 North in units of 1km and so on for increasing numbers of digits.

Not all grid squares contain land and as there is a false zero, the easiest way to picture it is as an image.

It is also possible to avoid all these complicated letters and numbers and simply quote two numbers, the distance East and distance North of the reference point, which is at the bottom left of the map. The distance is usually expressed in metres.

const double FIELD_WIDTH = 20.0;
const double FIELD_HEIGHT = 10.0;
const double SQUARE_WIDTH = 2.0;
const double SQUARE_HEIGHT = 1.0;
const double SUB_WIDTH = 0.0833333333333333;
const double SUB_HEIGHT = 0.0416666666666667;
const double SUB_SUB_WIDTH = 0.0083333333333333;
const double SUB_SUB_HEIGHT = 0.0041666666666667;

AnsiString DegToMaidenhead(double lat, double lon) // Returns Maidenhead Locator

{
char buff[9];

AnsiString locator;
double latitude = lat;
double longitude = lon;
sprintf(buff,"AA00AA00");

latitude += 90.0;
longitude += 180.0;

char ilong = char(longitude / FIELD_WIDTH);
buff[0] += ilong;
longitude -= double(ilong) * FIELD_WIDTH;

char ilat = char(latitude / FIELD_HEIGHT);
buff[1] += ilat;
latitude -= double(ilat) * FIELD_HEIGHT;

ilong = char(longitude / SQUARE_WIDTH);
buff[2] += ilong;
longitude -= double(ilong) * SQUARE_WIDTH;

ilat = char(latitude / SQUARE_HEIGHT);
buff[3] += ilat;
latitude -= double(ilat) * SQUARE_HEIGHT;

ilong = char(longitude / SUB_WIDTH);
buff[4] += ilong;
longitude -= double(ilong) * SUB_WIDTH;

ilat = char(latitude / SUB_HEIGHT);
buff[5] += ilat;
latitude -= double(ilat) * SUB_HEIGHT;

ilong = char(longitude / SUB_SUB_WIDTH);
buff[6] += ilong;

ilat = char(latitude / SUB_SUB_HEIGHT);
buff[7] += ilat;

locator.printf("%s",buff);
return locator;
}

MaidenheadToDeg(AnsiString locator, double &latitude, double &longitude)
{
// Converts Maidenhead to Lat Long, returns 1 on success

locator = locator.UpperCase();

if (locator.Length() == 4)
locator += "IL";

if (locator.Length() != 6 && locator.Length() != 8)
return 0;

if (locator[1] < 'A' || locator[1] > 'R' ||
locator[2] < 'A' || locator[2] > 'R' ||
locator[3] < '0' || locator[3] > '9' ||
locator[4] < '0' || locator[4] > '9' ||
locator[5] < 'A' || locator[5] > 'X' ||
locator[6] < 'A' || locator[6] > 'X')
return 0;

if (locator.Length() == 8 &&
(locator[7] <'0' || locator[7] > '9' ||
locator[8] < '0' || locator[8] > '9'))
return 0;

longitude = -180.0 + FIELD_WIDTH * (locator[1] - 'A') +
SQUARE_WIDTH * (locator[3] - '0') +
SUB_WIDTH * (locator[5] - 'A');

latitude = -90.0 + FIELD_HEIGHT * (locator[2] - 'A') +
SQUARE_HEIGHT * (locator[4] - '0') +
SUB_HEIGHT * (locator[6] - 'A');

if (locator.Length() == 8) {
longitude += SUB_SUB_WIDTH * (locator[7] - '0') + SUB_SUB_WIDTH / 2.0;
latitude += SUB_SUB_HEIGHT * (locator[8] - '0') + SUB_SUB_HEIGHT / 2.0;
} else {
longitude += SUB_SUB_WIDTH * 5.0;
latitude += SUB_SUB_HEIGHT * 5.0;
locator = locator + "55";
}
return 1;
}

const double deg2rad = 0.0174532925199432957692369076848;
const double rad2deg = 57.295779513082320876798154814105;

void OSGBtoLL(const double OSGBNorthing, const double OSGBEasting, const char* OSGBZone, double& Lat, double& Long )
{
//converts OSGB coords 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

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 RefEasting=0, RefNorthing=0;

//only calculate M0 once since it is based on the origin of the OSGB projection, which is fixed

static 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));

OSGBSquareToRefCoords(OSGBZone, RefEasting, RefNorthing);

x = OSGBEasting - 400000.0 + RefEasting; //remove 400,000 meter false easting for longitude
y = OSGBNorthing + 100000.0 + RefNorthing; //remove 100,000 meter false northing for longitude

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);

Lat = 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);

Lat = Lat * rad2deg;
Long = (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);
Long = LongOrigin + Long * rad2deg;
}

int LLtoOSGB(const double Lat, const double Long, long &OSGBEasting, long &OSGBNorthing, char OSGBGridSquare[3])
{
//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
// returns 1 if within UK, otherwise zero
//Written by Chuck Gantz- chuck.gantz@globalstar.com

// Fast sanity check - are we anywhere near the UK
if(Lat < 45 || Lat > 65)
return 0;
if(Long < -15 || Long > 5)
return 0;

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 = Lat*deg2rad;
double LongRad = Long*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
static 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

// Sanity check
if(easting < 0.0 || easting > 700000.0 || northing < 0.0 || northing > 1300000.0)
return(0);

CoordsToOSGBSquare(easting, northing, OSGBGridSquare, OSGBEasting, OSGBNorthing);

return(1); // Valid location
}

void OSGBSquareToRefCoords(const char* OSGBGridSquare, long &RefEasting, long &RefNorthing)
{
int pos, x_multiplier, y_multiplier;
char GridSquare[ ] = "VWXYZQRSTULMNOPFGHJKABCDE";

//find 500km offset
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;
}
RefEasting = x_multiplier * 500000L;
RefNorthing = y_multiplier * 500000L;

//find 100km offset and add to 500km offset to get coordinate of
//square point is in
pos = strchr(GridSquare, OSGBGridSquare[1]) - GridSquare;
RefEasting += ((pos % 5) * 100000L);
RefNorthing += ((pos / 5) * 100000L);

}

void CoordsToOSGBSquare(double easting, double northing, char OSGBGridSquare[3], long &OSGBEasting, long &OSGBNorthing)
{
char GridSquare[] = "VWXYZQRSTULMNOPFGHJKABCDE";
int posx, posy; //positions in grid

OSGBEasting = long(easting + 0.5); //round to nearest int
OSGBNorthing = long(northing + 0.5); //round to nearest int

//find correct 500km square
posx = OSGBEasting / 500000L;
posy = OSGBNorthing / 500000L;
OSGBGridSquare[0] = GridSquare[posx + posy * 5 + 7];

//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
OSGBGridSquare[1] = GridSquare[posx + posy * 5];

OSGBGridSquare[2] = '\0';//terminate grid ref string with null

//remainder is northing and easting
OSGBNorthing = OSGBNorthing % 500000L;
OSGBNorthing = OSGBNorthing % 100000L;

OSGBEasting = OSGBEasting % 500000L;
OSGBEasting = OSGBEasting % 100000L;
}