从EPSG 3857转换坐标到4326

9

我在数据库中有一组以EPSG 3857格式存储的坐标列表。 我需要将它们转换为EPSG 4326格式。 我尝试使用DotSpatial,但我的代码总是返回一个无限大的双精度数组。

public double[] ConvertCoodinates()
    {
        double[] xy = new double[2];
        xy[0] = 5085240.8300000000;
        xy[1] = 1530088.9600000000;
    //An array for the z coordinate
        double[] z = new double[1];
        z[0] = 0;
        ProjectionInfo pStart = KnownCoordinateSystems.Geographic.World.WGS1984;
        pStart.AuthorityCode = 3857;
        ProjectionInfo pEnd = KnownCoordinateSystems.Geographic.World.WGS1984;
        pEnd.AuthorityCode = 4326;
        Reproject.ReprojectPoints(xy, z, pStart, pEnd, 0, 1);
        return xy;
    }

这个xy数组总是包含无限大; 有人能帮忙吗?

4个回答

11

最后我找到了一个数学公式来转换这些坐标。

我将其实现为一个存储过程,因为我有一系列的点,而这个存储过程可以计算它们之间的距离。

DECLARE @e FLOAT=2.7182818284
DECLARE @X DECIMAL(18,2) =20037508.34

SET @StartLat3857 =(SELECT TOP 1 Latitude FROM Coordinates WHERE IdCoord=@IdCoord ORDER By IdTDFPath ASC)
SET @StartLng3857=(SELECT TOP 1 Longitude FROM Coordinates WHERE IdCoord=@IdCoord ORDER By IdTDFPath ASC)

--converting the logitute from epsg 3857 to 4326
            SET @StartLng=(@StartLng3857*180)/@X

--converting the latitude from epsg 3857 to 4326
            SET @StartLat = @StartLat3857/(@X/180)
            SET @StartLat = ((ATAN(POWER(@e,((PI()/180)*@StartLat))))/(PI()/360))-90

最终这只是一种数学公式,可以在任何语言中使用。例如,在JavaScript中它将是:

(a**2 + b**2) ** 0.5;

const e = 2.7182818284
const X = 20037508.34

const lat3857 = 1743704.947843 
const long3857 = 16978473.105100

//converting the logitute from epsg 3857 to 4326
const long4326 = (lat3857*180)/X

//converting the latitude from epsg 3857 to 4326 split in multiple lines for readability

let lat4326 = lat3857/(X / 180)
const exponent = (Math.PI / 180) * lat4326

lat4326 = Math.atan(e ** exponent)
lat4326 = lat4326 / (Math.PI / 360)
lat4326 = lat4326 - 90

2
你能帮我把4326转换回3857吗? - ABH
2
@ABH 你需要使用反向公式来计算坐标。 在C#中,它将是 private double[] ConvertCoordinate(double lat, double lng) { double x = lng * 20037508.34 / 180; double y = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180); y = y * 20037508.34 / 180; return new double[] { x, y }; } - Sethlans
@sethlans,你能否更新一下你的帖子。并不是每个人都理解存储过程。 - chitgoks
1
@chitgoks 希望 JavaScript 版本更易读 - Sethlans
我会分享我的Java版本。我刚刚做了一个。 - chitgoks
显示剩余3条评论

8

似乎选中答案的Javascript版本出现了错误。它总是给出比实际应该多两倍的答案,所以最终我选择将结果除以360而不是除以180(但我承认这里的数学超出了我的专业领域)。

以下是可行的代码:

  coord3857To4326(coord) {
    
    const e_value = 2.7182818284;
    const X = 20037508.34;
    
    const lat3857 = coord.lat
    const long3857 = coord.lng;
    
    //converting the longitute from epsg 3857 to 4326
    const long4326 = (long3857*180)/X;
    
    //converting the latitude from epsg 3857 to 4326 split in multiple lines for readability        
    let lat4326 = lat3857/(X / 180);
    const exponent = (Math.PI / 180) * lat4326;
    
    lat4326 = Math.atan(Math.pow(e_value, exponent));
    lat4326 = lat4326 / (Math.PI / 360); // Here is the fixed line
    lat4326 = lat4326 - 90;

    return {lat:lat4326, lng:long4326};
    
}

我也成功地进行了反向操作:

  coord4326To3857(coord) {

    const X = 20037508.34;

    let long3857 = (coord.lng * X) / 180;

    let lat3857 = parseFloat(coord.lat) + 90;
    lat3857 = lat3857 * (Math.PI/360);
    lat3857 = Math.tan(lat3857);
    lat3857 = Math.log(lat3857);
    lat3857 = lat3857 / (Math.PI / 180);
    
    lat3857 = (lat3857 * X) / 180;

    return {lat:lat3857, lng:long3857};
}

希望这能帮助到其他人!


2
谢谢,有一个错误,我已经在主要答案中进行了更正。 - Sethlans

4

EPSG 3857 坐标投影到 EPSG 4326 坐标系有些棘手。微软建议使用 ProjNet4GeoAPI,因此我决定使用它。

https://learn.microsoft.com/en-us/ef/core/modeling/spatial#srid-ignored-during-client-operations

我已验证在这里可以工作:

http://epsg.io/transform#s_srs=3857&t_srs=4326&x=1530088.9600000&y=5085240.8300000

示例转换:

var x = 1530088.96d;
var y = 5085240.83d;

var epsg3857ProjectedCoordinateSystem = ProjNet.CoordinateSystems.ProjectedCoordinateSystem.WebMercator;
var epsg4326GeographicCoordinateSystem = ProjNet.CoordinateSystems.GeographicCoordinateSystem.WGS84;

var coordinateTransformationFactory = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory();
var coordinateTransformation = coordinateTransformationFactory.CreateFromCoordinateSystems(epsg3857ProjectedCoordinateSystem, epsg4326GeographicCoordinateSystem);

var epsg3857Coordinate = new GeoAPI.Geometries.Coordinate(x, y);

var epsg4326Coordinate = coordinateTransformation.MathTransform.Transform(epsg3857Coordinate);

0

我使用以下代码在Rust中使其工作:

use std::f64::consts::{E, PI};

fn convert_coordinates(pos: &Vec<f64>) -> Vec<f64> {
    
    // Longitude <-> X
    // Latitude <-> Y
    // Output Decimal degree (WGS84)
    // X:-7920923.8631988168, Y: 5197114.9130207254
    //            0                     1
    // Longitude: -71.1548697°, Latitude: 42.2407752°
    //
    
    let long3857: f64 = pos[0];
    let lat3857: f64 = pos[1];

    let long4326: f64 = (long3857 * 180.0) / 20037508.34;
    let lat4326: f64 = (lat3857 * 180.0) / 20037508.34;

    let lat_rad: f64 = (lat4326 / 180.0) * PI;
    let lat4326_final: f64 = (2.0 * lat_rad.exp().atan() - PI / 2.0).to_degrees();

    vec![long4326, lat4326_final]
}

assert_eq!(-71.1548697_f64.round(), convert_coordinates(&vec![-7920923.8631988168, 5197114.9130207254])[0].round());
println!("Longitude Passed");
assert_eq!(42.2407752_f64.round(), convert_coordinates(&vec![-7920923.8631988168, 5197114.9130207254])[1].round());
println!("Latitude Passed");

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接