3

NEW ZEALAND MAP GRID座標 ieNorthingEastingWGS84座標 ieLatitudeとに変換したいLongitude

私はインターネットで検索しましたが、これを行う方法の適切な説明、またはこれを行うためのオンライン計算機はありません.

私の最終的な目標は、座標を座標に変換するプログラムをC#orで書くことです。JAVANZMGWGS84

4

3 に答える 3

2

注意: この実装では、ピンポイントの精度で変換することはできません。

これは、NZMG から NZGD1949 に変換し、次に NZGD1949 から NZGD2000 にデータム シフトするためにまとめた ac# 実装です。NZGD2000 座標は明らかに WGS84 ソースと互換性があります: http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/geodetic-datum-conversions/wgs84-nzgd2000

結果が一貫して緯度で約 40 メートルずれていることがわかったため、データム シフト メソッドに「useStaticCorrection」パラメータを追加しました。これは入力座標によって異なります。これはロジックの間違いによるものかもしれませんが、私が知る限り、すべて正しいようです (参照されているソースによると)。


出典:
定数と係数
http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/projection-conversions/new-zealand-map-grid

データム情報。
http://www.linz.govt.nz/regulatory/25000

NZMG -> NZGD1949 式
http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/projection-conversions/new-zealand-map-grid

NZ1949 -> NZGD2000 の違い
http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/geodetic-datum-conversions/nzgd1949-nzgd2000/three-seven

3 パラメーター シフトと同等です。
http://sas2.elte.hu/tg/eesti_datums_egs9.htm

//緯度/経度のファッジ量
http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/geodetic-datum-conversions/datum-transformation-examples


public class Pair<T, U>
{
    public Pair()
    {
    }

    public Pair(T first, U second)
    {
        this.First = first;
        this.Second = second;
    }

    public T First { get; set; }
    public U Second { get; set; }
};



using System.Numerics;

public class NZMGConverter
{
    //Constants and coefficients
    //Source: http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/projection-conversions/new-zealand-map-grid
    private static Double N0 = 6023150;
    private static Double E0 = 2510000;
    private static Double NZGD1949a = 6378388;
    private static Double NZMGLatOrigin = -41;
    private static Double NZMGLongOrigin = 173;
    //some of these are unused but are included for completeness. (unused arrays are for converting from wgs84 / nzgd2000 -> nzmg).
    private static Double[] A = new Double[] { 0.6399175073, -0.1358797613, 0.0632944090, -0.0252685300, 0.0117879000, -0.0055161000, 0.0026906000, -0.0013330000, 0.0006700000, -0.0003400000 };
    private static Complex[] B = new Complex[] { new Complex(0.7557853228, 0.0), new Complex(0.249204646, 0.003371507), new Complex(-0.001541739, 0.04105856), new Complex(-0.10162907, 0.01727609), new Complex(-0.26623489, -0.36249218), new Complex(-0.6870983, -1.1651967) };
    private static Complex[] C = new Complex[] { new Complex(1.3231270439, 0.0), new Complex(-0.577245789, -0.007809598), new Complex(0.508307513, -0.112208952), new Complex(-0.15094762, 0.18200602), new Complex(1.01418179, 1.64497696), new Complex(1.9660549, 2.5127645) };
    private static Double[] D = new Double[] { 1.5627014243, 0.5185406398, -0.0333309800, -0.1052906000, -0.0368594000, 0.0073170000, 0.0122000000, 0.0039400000, -0.0013000000 };

    //Datum info.
    //Source: http://www.linz.govt.nz/regulatory/25000
    private static Double NZGD1949f = 0.003367003;
    //Double NZGD1949InverseFlattening = 297;
    private static Double NZGD1949e2 = (2 * 0.003367003) - Math.Pow(0.003367003, 2);

    private static Double NZGD2000a = 6378137;
    private static Double NZGD2000f = 0.003352811;
    //Double NZGD2000InverseFlattening = 298.2572221;
    private static Double NZGD2000e2 = (2 * 0.003352811) - Math.Pow(0.003352811, 2);

    //NZ1949 -> NZGD2000 differences
    //Source: http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/geodetic-datum-conversions/nzgd1949-nzgd2000/three-seven
    private static Double differenceX = 54.4;
    private static Double differenceY = -20.1;
    private static Double differenceZ = 183.1;

    //Source: http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/projection-conversions/new-zealand-map-grid
    public static Pair<Double, Double> convertToNZGD1949(Double nzmgX, Double nzmgY)
    {
        Pair<Double, Double> result = null;

        Complex z = new Complex((nzmgY - N0) / NZGD1949a, (nzmgX - E0) / NZGD1949a);

        Complex theta0 = new Complex();
        for (int i = 0; i < C.Length; ++i) { theta0 += Complex.Multiply(C[i], Complex.Pow(z, i + 1)); }

        Double deltaLatWtheta0 = 0;
        for (int i = 0; i < D.Length; ++i) { deltaLatWtheta0 += D[i] * Math.Pow(theta0.Real, (i + 1)); }

        result = new Pair<double, double>();
        result.First = NZMGLatOrigin + (Math.Pow(10, 5) / 3600) * deltaLatWtheta0; //NZGD1949 lat
        result.Second = NZMGLongOrigin + 180 / Math.PI * theta0.Imaginary; //NZGD1949 long

        return result;
    }
    //Equivilent of a Three parameter shift.
    //Source: http://sas2.elte.hu/tg/eesti_datums_egs9.htm
    public static Pair<Double, Double> datumShiftNZGD1949toNZGD2000(Double nzgd1949Lat, Double nzgd1949Long, Boolean useStaticCorrection)
    {
        Pair<Double, Double> result = null;

        Double latInRaidans = nzgd1949Lat * (Math.PI / 180);
        Double lngInRaidans = nzgd1949Long * (Math.PI / 180);

        Double n = NZGD1949a / Math.Sqrt((1 - NZGD1949e2 * Math.Pow(Math.Sin(latInRaidans), 2)) * (3 / 2));
        Double m = NZGD1949a * ((1 - NZGD1949e2) / (1 - NZGD1949e2 * Math.Pow(Math.Sin(latInRaidans), 2) * (3 / 2)));

        Double difF = NZGD2000f - NZGD1949f;
        Double difA = NZGD2000a - NZGD1949a;

        Double changeInLat = (-1 * differenceX * Math.Sin(latInRaidans) * Math.Cos(lngInRaidans)
            - differenceY * Math.Sin(latInRaidans) * Math.Sin(lngInRaidans)
            + differenceZ * Math.Cos(latInRaidans) + (NZGD1949a * difF + NZGD1949f * difA) * Math.Sin(2 * latInRaidans))
            / (m * Math.Sin(1));
        Double changeInLong = (-1 * differenceX * Math.Sin(lngInRaidans) + differenceY * Math.Cos(lngInRaidans)) / (n * Math.Cos(latInRaidans) * Math.Sin(1));

        result = new Pair<Double, Double>();
        result.First = (latInRaidans + changeInLat) * (180 / Math.PI);
        result.Second = (lngInRaidans + changeInLong) * (180 / Math.PI);

        if (useStaticCorrection)
        {
            //Difference found by 3-parameter test values @ http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/geodetic-datum-conversions/datum-transformation-examples
            Double staticLatCorrection = 0.00037432;
            Double staticLongCorrection = -0.00003213;
            result.First -= staticLatCorrection;
            result.Second += staticLongCorrection;
        }

        return result;
    }
}
于 2014-12-05T00:14:50.003 に答える
1

Swiftを使用してNZMGから緯度/経度に。上記の DavidG の投稿に感謝します。緯度40mの誤差を修正。

import Foundation
import Darwin

class NZMGConverter {

    private let easting: Double!
    private let northing: Double!
    private var z: Complex! = nil
    private let N0: Double  = 6023150; //false Northing
    private let E0: Double  = 2510000; //false Easting

    private let NZMGLatOrigin: Double  = -41;
    private let NZMGLongOrigin: Double  = 173;

    private let c: [Complex]  = [Complex(real: 1.3231270439, imag: 0.0), Complex(real: -0.577245789, imag: -0.007809598), Complex(real: 0.508307513, imag: -0.112208952), Complex(real: -0.15094762, imag: 0.18200602), Complex(real: 1.01418179, imag: 1.64497696), Complex(real: 1.9660549, imag: 2.5127645)]
    private let b: [Complex] = [Complex(real: 0.7557853228, imag: 0.0), Complex(real: 0.249204646, imag: 0.003371507), Complex(real: -0.001541739, imag: 0.04105856), Complex(real: -0.10162907, imag: 0.01727609), Complex(real: -0.26623489, imag: -0.36249218), Complex(real: -0.6870983, imag: -1.1651967)]
    private let d: [Double] = [1.5627014243, 0.5185406398, -0.0333309800, -0.1052906000, -0.0368594000, 0.0073170000, 0.0122000000, 0.0039400000, -0.0013000000]

    //Datum info.
    private let NZGD1949f: Double  = 0.003367003;
    //Double NZGD1949InverseFlattening = 297;
    private let NZGD1949e2: Double = (2 * 0.003367003) - pow(0.003367003, 2);
    private let NZGD1949a: Double  = 6378388; //Semi-major axis

    private let NZGD2000a: Double  = 6378137;
    private let NZGD2000f: Double  = 0.003352811;
    //Double NZGD2000InverseFlattening = 298.2572221;
    private let NZGD2000e2: Double  = (2 * 0.003352811) - pow(0.003352811, 2);

    //NZ1949 -> NZGD2000 differences
    private let differenceX: Double = 54.4;
    private let differenceY: Double  = -20.1;
    private let differenceZ: Double = 183.1;

    init(easting: Double, northing: Double) {

        self.easting = easting
        self.northing = northing
    }

    func nzmgToNZGD1949() -> (latitude: Double, longitude: Double) {

        z = Complex(real: (northing - N0) / NZGD1949a, imag: (easting - E0) / NZGD1949a)
        let theta0 = thetaZero(z)
        let theta: Complex = thetaSuccessive(theta0, i: 0)
        let latlong = calcLatLong(theta)
        print(latlong)
        return datumShift(latlong.0, longitude: latlong.1)
    }

    func thetaZero(z: Complex) -> Complex {
        var theta0 = Complex()
        for var i = 0; i < c.count; i++ {
            theta0 = theta0 + (c[i] * (z ^ (i + 1)))
            //print(theta0)
        }
        return theta0
    }

    func thetaSuccessive(theta0: Complex, i: Int) -> Complex {
        var numerator: Complex = Complex()
        var denominator: Complex = Complex()
        for var n = 1; n < b.count; n++ {
            numerator = numerator + (Complex(real: n) * b[n] * (theta0 ^ (n + 1)))
        }
        numerator = z + numerator
        for var n = 0; n < b.count; n++ {
            denominator = denominator + (Complex(real: n + 1) * b[n] * (theta0 ^ n))
        }
        let theta = numerator / denominator
        if i == 2 {
            return theta
        } else {
            return thetaSuccessive(theta, i: i + 1)
        }
    }

    func calcLatLong(theta: Complex) -> (latitude: Double, longitude: Double) {

        let deltaPsi = theta.real
        let deltaLambda = theta.imag
        var deltaPhi: Double = 0
        for var n = 0; n < d.count; n++ {
            deltaPhi += d[n] * pow(deltaPsi, Double(n + 1))
        }
        let lat: Double = NZMGLatOrigin + ((pow(10, 5) / 3600) * deltaPhi)
        let long: Double = NZMGLongOrigin + ((180 / M_PI) * deltaLambda)
        return (latitude: lat, longitude: long)
    }

    func datumShift(latitude: Double, longitude: Double) -> (latitude: Double, longitude: Double) {

        let latInRadians = latitude * (M_PI / 180)
        let lngInRadians = longitude * (M_PI / 180)

        let v = NZGD1949a / (sqrt(1 - (NZGD1949e2 * pow(sin(latInRadians), 2))))
        var x_cart = v * cos(latInRadians) * cos(lngInRadians)
        var y_cart = v * cos(latInRadians) * sin(lngInRadians)
        var z_cart = (v * (1 - NZGD1949e2)) * sin(latInRadians)
        x_cart += differenceX
        y_cart += differenceY
        z_cart += differenceZ
        print(x_cart)
        print(y_cart)
        print(z_cart)

        let p = sqrt(pow(x_cart, 2) + pow(y_cart, 2))
        let r = sqrt(pow(p, 2) + pow(z_cart, 2))
        let mu = atan((z_cart / p) * ((1 - NZGD2000f) + ((NZGD2000e2 * NZGD2000a) / r)))

        let num = (z_cart * (1 - NZGD2000f)) + (NZGD2000e2 * NZGD2000a * pow(sin(mu), 3))
        let denom = (1 - NZGD2000f) * (p - (NZGD2000e2 * NZGD2000a * pow(cos(mu), 3)))

        var lat = atan(num / denom)
        var long = atan(y_cart / x_cart)
        print(long)

        lat = (180 / M_PI) * lat
        long = 180 + (180 / M_PI) * long

        return (latitude: lat, longitude: long)

    }


}
}
于 2015-10-30T01:03:51.880 に答える
0

Proj4j は次の機能を提供します。これは、数少ないユニバーサル座標変換ライブラリの 1 つです。 http://trac.osgeo.org/proj4j/

次に、新しい熱意とデータムの説明を定義する proj4 文字列を見つける必要があります。

古いグリッドが必要か新しいグリッドが必要かによって、proj4 に含まれているニュージーランドのデータム シフト グリッドが必要になります。

最初のステップとして、必要なデータムの名前を正確に見つけます: NZGD49 または NZGD2000 ?

于 2013-08-14T15:57:54.047 に答える