是否有任何已知的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;
}
看看这页
在"程序/源代码"一节中,您将发现题为"WGS84 Lat/Long <=> OSGB36 Grid References"的小节。在本小节中,有几个到提供此功能的实用程序的链接。
在另一篇文章中,你可以找到这样一个任务的源代码和解释。
如果您不确定其他库的正确性,我会将结果与PROJ.4进行比较。如果你不想使用本机库,也有一个纯Java的PROJ.4端口。
使用project .4用于WGS84 l<=> OSGB36已经在这里讨论过了:
project .4库和OSGB36
:
http://trac.osgeo.org/proj/wiki/man_cs2cs
OSGB36字符串:
http://spatialreference.org/ref/epsg/27700/proj4/
WGS84字符串:
http://spatialreference.org/ref/epsg/4326/proj4/
Java端口:http://sourceforge.net/projects/jmapprojlib/
我在这方面做了一些工作,并使用了文档《英国坐标系统指南》,如果你仔细阅读,它会告诉你所有你需要知道的。从公式中使用的变量名来看,您很可能已经研究了同一文档。我认为整个文档中关于坐标系转换最重要的段落是这段(从第30页底部开始)
总结:对于一个简单的基准,纬度和经度坐标从基准a到基准B, 首先将转换为笛卡尔坐标(公式见附件B),取所有椭球体高度为零点,利用基准面A的椭球参数;则从基准a应用Helmert变换用式(3)对基准B;finally使用椭球体转换回经纬度基准B的参数(公式见附件C),丢弃基准B椭球高度。
它描述了你必须采取的三个步骤。您引用的公式将把纬度/经度转换为同一基准中的东/北。它不是一个变换
从你的方法的命名中,你正在传递WGS84的纬度/经度,所以你应该:
1)把网格参考系统(以及相关的真实原点等)的所有想法都抛在脑后,直到你把经纬度从一个基准转换到另一个基准为止
2)使用操作系统指南中的公式B1到B5,将WGS84的经纬度转换为WGS84基准的3D笛卡尔坐标(即x, y和z)。确保为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基准面的经纬度。
从这里开始,你可以使用你引用的公式来计算网格的东北方
我使用jcoord。