C++ Convert Lat Long to BNG with Proj.4

Question

Im creating a plugin for VBS2. In the program I can output a grid in either LAT LONG,UTM or MGRS. I need to be able to convert to BNG. I have managed to create a TKinter application in python that works using Proj.4 but now need to create it in C++ as a DLL.

The LatLong im using (51.20650N 1.81906W) is a known point, the BNG conversion is approx SU 127 452.

#include <proj_api.h>




double x = 51.20650; //atof(GetStrArgument(input, 0)); 
double y = 1.81906;  //atof(GetStrArgument(input, 1));  

char *pj_latlongc = "+proj=longlat +datum=WGS84";
char *pj_UTMc = "+proj=utm +zone=29 +ellps=WGS84";
char *osc = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs";

projPJ pj_OS = pj_init_plus(osc);
projPJ pj_UTM = pj_init_plus(pj_UTMc);
projPJ pj_latlong = pj_init_plus(pj_latlongc);

x *= DEG_TO_RAD;
y *= DEG_TO_RAD;

int p = pj_transform(pj_latlong, pj_OS, 1, 1, &x, &y, NULL);

Unfortunately the result is completely wrong and doesn't make sense.Could someone help on this?

Answer 1

Switching x for y and flipping the sign of x seems to yield better results. Running those numbers on the console, the proj output seems to agree better with the point you specified (51.20650N 1.81906W)

echo -1.81906 51.20650 | proj -V +proj=tmerc
Longitude: 1d49'8.616"W [ -1.81906 ]
Latitude:  51d12'23.4"N [ 51.2065 ]

Similarly, cs2cs returns seemingly sensible output:

echo -1.81906 51.20650 | cs2cs -v +proj=longlat +datum=WGS84   +to   +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs
# ---- From Coordinate System ----
#Lat/long (Geodetic alias)
#
# +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# ---- To Coordinate System ----
#Transverse Mercator
# Cyl, Sph&Ell
# +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000
# +ellps=airy +datum=OSGB36 +units=m +no_defs
# +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894
#--- following specified but NOT used
# +ellps=airy
412736.92  145270.05 -47.95

This is identical to the C++ program's output

double x = -1.81906;
double y = 51.20650;
412736.924409  145270.054358

For comparison, I ran the numbers through the BNG converter of the British Geological Survey, which seems to agree as well

Easting: 412737
Northing: 145270

For writing these numbers in the form of an Ordnance Survey National Grid reference, just subtract integer multiples of 100km to derive the corresponding letter as illustrated in this Source Code for JavaScript. The final result is SU 127 452, as desired:

#include <proj_api.h>
#include <iostream>

int main() {
  double x = -1.81906;
  double y = 51.20650;

  const char* pj_latlongc = "+proj=longlat +datum=WGS84";
  const char* osc = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs";

  projPJ pj_latlong = pj_init_plus(pj_latlongc);
  projPJ pj_OS = pj_init_plus(osc);

  x *= DEG_TO_RAD;
  y *= DEG_TO_RAD;

  int p = pj_transform(pj_latlong, pj_OS, 1, 1, &x, &y, NULL);

  std::cout.setf(std::ios::fixed, std::ios::floatfield);
  std::cout.setf(std::ios::showpoint);
  std::cout << x << "  " << y << std::endl;;
  // prints 412736.924409  145270.054358


  // now convert to UK Grid Ref
  int e1 = floor(x/100000);
  int n1 = floor(y/100000);
  int e2 = (int)x % 100000 / 100;
  int n2 = (int)y % 100000 / 100;
  char l1 = (19 - n1) - (19-n1) % 5 + ((e1 + 10)/5);
  char l2 = (19 - n1) * 5 % 25 + e1 % 5;
  if (l1 > 7) l1++;
  if (l2 > 7) l2++;
  l1 = 'A' + l1;
  l2 = 'A' + l2;

  std::cout << l1 << l2 << ' ' << e2 << ' ' << n2 << std::endl;
  // prints SU 127 452
}