經緯度轉換TWD97

from math import tan, sin, cos, radians
 
class LatLonToTWD97(object):
    """This object provide method for converting lat/lon coordinate to TWD97
    coordinate
 
    the formula reference to
    http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm (there is lots of typo)
    http://www.offshorediver.com/software/utm/Converting UTM to Latitude and Longitude.doc
 
    Parameters reference to
    http://rskl.geog.ntu.edu.tw/team/gis/doc/ArcGIS/WGS84%20and%20TM2.htm
    http://blog.minstrel.idv.tw/2004/06/taiwan-datum-parameter.html
    """
 
    def __init__(self,
        a = 6378137.0,
        b = 6356752.314245,
        long0 = radians(121),
        k0 = 0.9999,
        dx = 250000,
    ):
        # Equatorial radius
        self.a = a
        # Polar radius
        self.b = b
        # central meridian of zone
        self.long0 = long0
        # scale along long0
        self.k0 = k0
        # delta x in meter
        self.dx = dx
 
    def convert(self, lat, lon):
        """Convert lat lon to twd97 
 
        """
        a = self.a
        b = self.b
        long0 = self.long0
        k0 = self.k0
        dx = self.dx
 
        e = (1-b**2/a**2)**0.5
        e2 = e**2/(1-e**2)
        n = (a-b)/(a+b)
        nu = a/(1-(e**2)*(sin(lat)**2))**0.5
        p = lon-long0
 
        A = a*(1 - n + (5/4.0)*(n**2 - n**3) + (81/64.0)*(n**4  - n**5))
        B = (3*a*n/2.0)*(1 - n + (7/8.0)*(n**2 - n**3) + (55/64.0)*(n**4 - n**5))
        C = (15*a*(n**2)/16.0)*(1 - n + (3/4.0)*(n**2 - n**3))
        D = (35*a*(n**3)/48.0)*(1 - n + (11/16.0)*(n**2 - n**3))
        E = (315*a*(n**4)/51.0)*(1 - n)
 
        S = A*lat - B*sin(2*lat) + C*sin(4*lat) - D*sin(6*lat) + E*sin(8*lat)
 
        K1 = S*k0
        K2 = k0*nu*sin(2*lat)/4.0
        K3 = (k0*nu*sin(lat)*(cos(lat)**3)/24.0) * \
            (5 - tan(lat)**2 + 9*e2*(cos(lat)**2) + 4*(e2**2)*(cos(lat)**4))
 
        y = K1 + K2*(p**2) + K3*(p**4)
 
        K4 = k0*nu*cos(lat)
        K5 = (k0*nu*(cos(lat)**3)/6.0) * \
            (1 - tan(lat)**2 + e2*(cos(lat)**2))
 
        x = K4*p + K5*(p**3) + self.dx
        return x, y
if __name__ == '__main__':
    from math import degrees
 
    c = LatLonToTWD97()
    lat = radians(float(raw_input('lat:')))
    lon = radians(float(raw_input('lon:')))
    print 'input lat/lon', degrees(lat), degrees(lon)
    x, y = c.convert(lat, lon)
    print (x, y)

package com.ez2learn.android.powergrid.geo;
 
public class TMToLatLon {
	/**
		@brief TM coordinate to lat lon formula
 
		reference to http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm
 
		@param TMParameter object
		@param x
		@param y
	**/
	static public double[] convert(TMParameter tm, double x, double y) {
		double dx = tm.getDx();
		double dy = tm.getDy();
		double lon0 = tm.getLon0();
		double k0 = tm.getK0();
		double a = tm.getA();
		double b = tm.getB();
		double e = tm.getE();
 
		x -= dx;
		y -= dy;
 
		// Calculate the Meridional Arc
		double M = y/k0;
 
		// Calculate Footprint Latitude
		double mu = M/(a*(1.0 - Math.pow(e, 2)/4.0 - 3*Math.pow(e, 4)/64.0 - 5*Math.pow(e, 6)/256.0));
		double e1 = (1.0 - Math.pow((1.0 - Math.pow(e, 2)), 0.5)) / (1.0 + Math.pow((1.0 - Math.pow(e, 2)), 0.5));
 
		double J1 = (3*e1/2 - 27*Math.pow(e1, 3)/32.0);
		double J2 = (21*Math.pow(e1, 2)/16 - 55*Math.pow(e1, 4)/32.0);
		double J3 = (151*Math.pow(e1, 3)/96.0);
		double J4 = (1097*Math.pow(e1, 4)/512.0);
 
		double fp = mu + J1*Math.sin(2*mu) + J2*Math.sin(4*mu) + J3*Math.sin(6*mu) + J4*Math.sin(8*mu);
 
		// Calculate Latitude and Longitude
 
		double e2 = Math.pow((e*a/b), 2);
		double C1 = Math.pow(e2*Math.cos(fp), 2);
		double T1 = Math.pow(Math.tan(fp), 2);
		double R1 = a*(1-Math.pow(e, 2))/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), (3.0/2.0));
		double N1 = a/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), 0.5);
 
		double D = x/(N1*k0);
 
		// lat
		double Q1 = N1*Math.tan(fp)/R1;
		double Q2 = (Math.pow(D, 2)/2.0);
		double Q3 = (5 + 3*T1 + 10*C1 - 4*Math.pow(C1, 2) - 9*e2)*Math.pow(D, 4)/24.0;
		double Q4 = (61 + 90*T1 + 298*C1 + 45*Math.pow(T1, 2) - 3*Math.pow(C1, 2) - 252*e2)*Math.pow(D, 6)/720.0;
		double lat = fp - Q1*(Q2 - Q3 + Q4);
 
		// long
		double Q5 = D;
		double Q6 = (1 + 2*T1 + C1)*Math.pow(D, 3)/6;
		double Q7 = (5 - 2*C1 + 28*T1 - 3*Math.pow(C1, 2) + 8*e2 + 24*Math.pow(T1, 2))*Math.pow(D, 5)/120.0;
		double lon = lon0 + (Q5 - Q6 + Q7)/Math.cos(fp);
 
		return new double[] {Math.toDegrees(lat), Math.toDegrees(lon)};
	}
}