C++ 使用 Proj.4 将 Lat Long 转换为 BNG

Question

我正在为 VBS2 创建一个插件。在程序中，我可以输出 LAT LONG、UTM 或 MGRS 格式的网格。我需要能够转换为 BNG。我已经设法在 python 中创建了一个使用 Proj.4 工作的 TKinter 应用程序，但现在需要在 C++ 中将其创建为 DLL。

LatLong im 使用 (51.20650N 1.81906W) 是已知点，BNG 转换约为 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);

不幸的是，结果是完全错误的，没有意义。有人可以帮忙吗？

Answer 1

将 x 换成 y 并翻转 x 的符号似乎会产生更好的结果。在控制台上运行这些数字，proj 输出似乎更符合 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 ]

类似地，cs2cs 返回看似合理的输出:

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

这与 C++ 程序的输出相同

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

为了比较，我通过 BNG converter of the British Geological Survey 运行了数字, 这似乎也同意

Easting: 412737
Northing: 145270

将这些数字写成 Ordnance Survey National Grid 的形式引用，只需减去 100km 的整数倍即可得出相应的字母，如图 Source Code for JavaScript .最终结果是SU 127 452，如愿:

#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
}