将 WGS84 转换为 OSGB36

发布于 2024-12-21 16:17:12 字数 2749 浏览 2 评论 0原文

是否有任何已知的 java 库允许我将 WGS 84 线转换为 OSGB36 或者是否有一个好的公式我可以使用？我目前正在使用这个，但它不是很准确，所以我想知道是否有更好的可以使用。

private static double[] Wgs84ToBNG(double inLat, double inLon) {

    double lat = inLat * Math.PI / 180.0;
    double lon = inLon * Math.PI / 180.0;
    double a = 6377563.396; // Airy 1830 major & minor semi-axes
    double b = 6356256.910;
    double F0 = 0.9996012717; // NatGrid scale factor on central meridian
    double lat0 = 49 * Math.PI / 180.0; // NatGrid true origin
    double lon0 = -2 * Math.PI / 180.0;
    double N0 = -100000; // northing & easting of true origin, metres
    double E0 = 400000;
    double e2 = 1 - (b * b) / (a * a); // eccentricity squared
    double n = (a - b) / (a + b), n2 = n * n, n3 = n * n * n;
    double cosLat = Math.cos(lat), sinLat = Math.sin(lat);
    double nu = a * F0 / Math.sqrt(1 - e2 * sinLat * sinLat); // transverse
                                                                // radius of
                                                                // curvature
    double rho = a * F0 * (1 - e2)
            / Math.pow(1 - e2 * sinLat * sinLat, 1.5); // meridional radius
                                                        // of curvature
    double eta2 = nu / rho - 1;
    double Ma = (1 + n + (5 / 4) * n2 + (5 / 4) * n3) * (lat - lat0);
    double Mb = (3 * n + 3 * n * n + (21 / 8) * n3) * Math.sin(lat - lat0)
            * Math.cos(lat + lat0);
    double Mc = ((15 / 8) * n2 + (15 / 8) * n3)
            * Math.sin(2 * (lat - lat0)) * Math.cos(2 * (lat + lat0));
    double Md = (35 / 24) * n3 * Math.sin(3 * (lat - lat0))
            * Math.cos(3 * (lat + lat0));
    double M = b * F0 * (Ma - Mb + Mc - Md); // meridional arc
    double cos3lat = cosLat * cosLat * cosLat;
    double cos5lat = cos3lat * cosLat * cosLat;
    double tan2lat = Math.tan(lat) * Math.tan(lat);
    double tan4lat = tan2lat * tan2lat;
    double I = M + N0;
    double II = (nu / 2) * sinLat * cosLat;
    double III = (nu / 24) * sinLat * cos3lat * (5 - tan2lat + 9 * eta2);
    double IIIA = (nu / 720) * sinLat * cos5lat
            * (61 - 58 * tan2lat + tan4lat);
    double IV = nu * cosLat;
    double V = (nu / 6) * cos3lat * (nu / rho - tan2lat);
    double VI = (nu / 120)
            * cos5lat
            * (5 - 18 * tan2lat + tan4lat + 14 * eta2 - 58 * tan2lat * eta2);
    double dLon = lon - lon0;
    double dLon2 = dLon * dLon;
    double dLon3 = dLon2 * dLon;
    double dLon4 = dLon3 * dLon;
    double dLon5 = dLon4 * dLon;
    double dLon6 = dLon5 * dLon;
    double N = I + II * dLon2 + III * dLon4 + IIIA * dLon6;
    double E = E0 + IV * dLon + V * dLon3 + VI * dLon5;
    double[] returnValue = { E, N };
    return returnValue;
}

原文

Is there any known java library which allows me to convert WGS 84 cords to OSGB36 or if there is a good formula I can use? I am currently using this one, but it's not very accurate so I'm wondering if there's a better one I can use.

private static double[] Wgs84ToBNG(double inLat, double inLon) {

    double lat = inLat * Math.PI / 180.0;
    double lon = inLon * Math.PI / 180.0;
    double a = 6377563.396; // Airy 1830 major & minor semi-axes
    double b = 6356256.910;
    double F0 = 0.9996012717; // NatGrid scale factor on central meridian
    double lat0 = 49 * Math.PI / 180.0; // NatGrid true origin
    double lon0 = -2 * Math.PI / 180.0;
    double N0 = -100000; // northing & easting of true origin, metres
    double E0 = 400000;
    double e2 = 1 - (b * b) / (a * a); // eccentricity squared
    double n = (a - b) / (a + b), n2 = n * n, n3 = n * n * n;
    double cosLat = Math.cos(lat), sinLat = Math.sin(lat);
    double nu = a * F0 / Math.sqrt(1 - e2 * sinLat * sinLat); // transverse
                                                                // radius of
                                                                // curvature
    double rho = a * F0 * (1 - e2)
            / Math.pow(1 - e2 * sinLat * sinLat, 1.5); // meridional radius
                                                        // of curvature
    double eta2 = nu / rho - 1;
    double Ma = (1 + n + (5 / 4) * n2 + (5 / 4) * n3) * (lat - lat0);
    double Mb = (3 * n + 3 * n * n + (21 / 8) * n3) * Math.sin(lat - lat0)
            * Math.cos(lat + lat0);
    double Mc = ((15 / 8) * n2 + (15 / 8) * n3)
            * Math.sin(2 * (lat - lat0)) * Math.cos(2 * (lat + lat0));
    double Md = (35 / 24) * n3 * Math.sin(3 * (lat - lat0))
            * Math.cos(3 * (lat + lat0));
    double M = b * F0 * (Ma - Mb + Mc - Md); // meridional arc
    double cos3lat = cosLat * cosLat * cosLat;
    double cos5lat = cos3lat * cosLat * cosLat;
    double tan2lat = Math.tan(lat) * Math.tan(lat);
    double tan4lat = tan2lat * tan2lat;
    double I = M + N0;
    double II = (nu / 2) * sinLat * cosLat;
    double III = (nu / 24) * sinLat * cos3lat * (5 - tan2lat + 9 * eta2);
    double IIIA = (nu / 720) * sinLat * cos5lat
            * (61 - 58 * tan2lat + tan4lat);
    double IV = nu * cosLat;
    double V = (nu / 6) * cos3lat * (nu / rho - tan2lat);
    double VI = (nu / 120)
            * cos5lat
            * (5 - 18 * tan2lat + tan4lat + 14 * eta2 - 58 * tan2lat * eta2);
    double dLon = lon - lon0;
    double dLon2 = dLon * dLon;
    double dLon3 = dLon2 * dLon;
    double dLon4 = dLon3 * dLon;
    double dLon5 = dLon4 * dLon;
    double dLon6 = dLon5 * dLon;
    double N = I + II * dLon2 + III * dLon4 + IIIA * dLon6;
    double E = E0 + IV * dLon + V * dLon3 + VI * dLon5;
    double[] returnValue = { E, N };
    return returnValue;
}

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

风追烟花雨 2024-12-28 16:17:12

看看这个页面

在“程序/源代码”部分，您会发现标题为“WGS84 纬度/经度 <=> OSGB36 网格参考”的小节。在本小节中，有几个指向提供此功能的实用程序的链接。

在另一篇帖子中，您可以找到以下内容的源代码这样的任务带有解释。

回复收藏 0 原文

落在眉间の轻吻 2024-12-28 16:17:12

如果您不确定其他库的正确性，我会将结果与 PROJ.4 进行比较。如果您不想使用本机库，还有一个用于 PROJ.4 的纯 Java 端口。

使用 PROJ.4 进行 WGS84 纬度/经度 <=> OSGB36 已在此讨论：
PROJ.4 库和 OSGB36

使用 cs2cs：
http://trac.osgeo.org/proj/wiki/man_cs2cs
OSGB36 字符串：
http://spatialreference.org/ref/epsg/27700/proj4/< br>
WGS84 字符串：
http://spatialreference.org/ref/epsg/4326/proj4/

Java港口：
http://sourceforge.net/projects/jmapprojlib/

回复收藏 0 原文

信愁 2024-12-28 16:17:12

我在这方面做了一些工作，并使用了文档指南到英国的坐标系如果你仔细阅读它，它会告诉你你需要知道的一切。
从公式中使用的变量名称来看，您很可能已经研究过同一文档。我认为整个文档中与坐标系之间转换相关的最重要的一段是（从第 30 页底部开始）

总结一下：对于纬度和经度坐标从基准 A 到
基准面B，首先转换为笛卡尔坐标（公式见附件B），取所有椭球体高度为
零并使用基准 A 的椭球参数；然后从数据 A 应用 Helmert 变换
使用等式(3)计算数据B；最后使用椭球体转换回纬度和经度
基准B的参数（附件C中的公式），丢弃基准B椭球高度。

它描述了您必须采取的3个步骤。您引用的公式会将纬度/经度转换为同一数据中的东距/北距。这不是一个转换

从您传递 WGS84 纬度/经度的方法的命名来看，您应该：

1）将网格参考系统（以及相关的真实原点等）的所有想法从您的脑海中剔除，直到您转换了纬度/lon 从一个基准到另一个基准

2) 使用公式 B1 将 WGS84 纬度/经度转换为 WGS84 基准的 3D 笛卡尔坐标（即 x、y 和 z）操作系统指南中的 B5。确保使用 WGS84 基准的参数（长轴/短轴）

3) 使用提到的 Helmert 变换，将刚刚计算的笛卡尔坐标转换为相对于 Airy 1830 椭球体的笛卡尔坐标。您将在第 6.6 节中找到获得新笛卡尔坐标所需的 7 个参数

新坐标 xGB、yGB、zGB 是：

double xGB = tx + (x * (1 + s)) + (-rz * y) + (ry * z);
double yGB = ty + (rz * x) + (y * (1 + s)) + (-rx * z);
double zGB = tz + (-ry * x) + (rx * y) + (z * (1 + s));

他们实际上并没有拼写出来，他们只是假设您记住了矩阵数学

4) 现在转换这些新坐标使用公式 B6 到 B8 将笛卡尔坐标（相对于 OSGB36 基准）转换为相对于 OSGB36 基准的纬度/经度

从这里您可以继续计算网格东距和使用您引用的公式进行北距

I've done a bit of work on this and have used the document A guide to coordinate systems in Great Britain It tells you all you need to know if you read it carefully.
From the look of the variable names used in your formula you may well have studied the same document. I think the most important paragraph in the whole document relating to converting between coordinate systems is this (from the bottom of page 30)

To summarise: For a simple datum change of latitude and longitude coordinates from datum A to
datum B, first convert to Cartesian coordinates (formulae in annexe B), taking all ellipsoid heights as
zero and using the ellipsoid parameters of datum A; then apply a Helmert transformation from datum A
to datum B using equation (3); finally convert back to latitude and longitude using the ellipsoid
parameters of datum B (formulae in annexe C), discarding the datum B ellipsoid height.

It describes the 3 steps you have to take. The formula you have quoted will convert a latitude/longitude to an easting/northing in the same datum. It's not a transformation

From the naming of your method you are passing WGS84 lat/lon in, so you should:

1) Put all thoughts of a grid reference system (and associated true origins etc) out of your mind until you have converted the lat/lon from one datum to the other

2) Convert that WGS84 lat/lon into the 3D Cartesians (that's the x, y and z) for the WGS84 datum, using the formulae B1 to B5 in the OS Guide. Make sure you use the parameters (major/minor axes) for the WGS84 datum

3) Using the Helmert transformation referred to, convert the Cartesians you've just calculated into Cartesians relative to the Airy 1830 ellipsoid. You'll find the 7 parameters you need to get the new Cartesians in section 6.6

The new coordinates xGB, yGB, zGB are:

double xGB = tx + (x * (1 + s)) + (-rz * y) + (ry * z);
double yGB = ty + (rz * x) + (y * (1 + s)) + (-rx * z);
double zGB = tz + (-ry * x) + (rx * y) + (z * (1 + s));

They don't actually spell that out, they just assume you remember your matrix maths

4) Now convert those new Cartesians (which are relative to the OSGB36 datum) into a lat/lon relative to that OSGB36 datum using formulae B6 to B8

From here you can go on to work out grid eastings and northings using that formula you have quoted

回复收藏 0 原文