Hej Engång i tiden för jättelänge sen hackade jag ihop följande för ett annat projekt som rann ut i sanden.
Borde gå att återanvända endel för att få till SWEREF99 13 30 00 package org.openstreetmap.josm.data.projection; import org.openstreetmap.josm.data.coor.EastNorth; import org.openstreetmap.josm.data.coor.LatLon; /** * Implemented from: * http://www.lantmateriet.se/upload/filer/kartor/geodesi_gps_och_detaljmatning/geodesi/Formelsamling/Gauss_Conformal_Projection.pdf */ public abstract class GaussKruger implements Projection { private Ellipsoid ellipsoid; // Central Meridian private double lambda_zero; // False Northing private double false_n; // False Easting private double false_e; // Scale factor private double k_zero; //Derived variables from the above private double f; private double n; private double a_roof; GaussKruger(Ellipsoid e, double central_meridian, double false_northing, double false_easting, double scale_factor) { ellipsoid = e; lambda_zero = Math.toRadians(central_meridian); false_n = false_northing; false_e = false_easting; k_zero = scale_factor; f = 1.0 - (ellipsoid.b / ellipsoid.a); n = f / (2.0 - f); a_roof = (ellipsoid.a / (1 + n)) * (1 + n * n / 4 + n * n * n * n / 64); a = ellipsoid.e2; b = (5 * a * a - Math.pow(a, 3)) / 6; c = (104 * Math.pow(a, 3) - 45 * Math.pow(a, 4)) / 120; d = 1237 * Math.pow(a, 4) / 1260; beta1 = n / 2 - 2 * n * n / 3 + 5 * Math.pow(n, 3) / 16 + 41 * Math.pow(n, 4) / 180; beta2 = 13 * n * n / 48 - 3 * Math.pow(n, 3) / 5 + 557 * Math.pow(n, 4) / 1440; beta3 = 61 * Math.pow(n, 3) / 240 - 103 * Math.pow(n, 4) / 140; beta4 = 49561 * Math.pow(n, 4) / 161280; astar = a + Math.pow(a, 2) + Math.pow(a, 3) + Math.pow(a, 4); bstar = (7 * Math.pow(a, 2) + 17 * Math.pow(a, 3) + 30 * Math.pow(a, 4)) / -6; cstar = (224 * Math.pow(a, 3) + 889 * Math.pow(a, 4)) / 120; dstar = 4279 * Math.pow(a, 4) / 1260; delta1 = n / 2 - 2 * n * n / 3 + 37 * Math.pow(n, 3) / 96 - Math.pow(n, 4) / 360; delta2 = n * n / 48 + Math.pow(n, 3) / 15 - 437 * Math.pow(n, 4) / 1440; delta3 = 17 * Math.pow(n, 3) / 480 - 37 * Math.pow(n, 4) / 840; delta4 = 4397 * Math.pow(n, 4) / 161280; } private GaussKruger() { } // latlon2eastNorth() variables private double a; private double b; private double c; private double d; private double beta1; private double beta2; private double beta3; private double beta4; //eastNorth2latlon() variables private double astar; private double bstar; private double cstar; private double dstar; private double delta1; private double delta2; private double delta3; private double delta4; public EastNorth latlon2eastNorth(LatLon p) { double lambda = Math.toRadians(p.lon()); double phi = Math.toRadians(p.lat()); double phi_star = phi - Math.sin(phi) * Math.cos(phi) * (a + b * Math.pow(Math.sin(phi), 2) + c * Math.pow(Math.sin(phi), 4) + d * Math.pow(Math.sin(phi), 6)); double delta_lambda = lambda - lambda_zero; double xi_prim = Math.atan(Math.tan(phi_star) / Math.cos(delta_lambda)); double eta_prim = atanh(Math.cos(phi_star) * Math.sin(delta_lambda)); double x = false_n + k_zero * a_roof * (xi_prim + beta1 * Math.sin(2 * xi_prim) * Math.cosh(2 * eta_prim) + beta2 * Math.sin(4 * xi_prim) * Math.cosh(4 * eta_prim) + beta3 * Math.sin(6 * xi_prim) * Math.cosh(6 * eta_prim) + beta4 * Math.sin(8 * xi_prim) * Math.cosh(8 * eta_prim)); double y = false_e + k_zero * a_roof * (eta_prim + beta1 * Math.cos(2 * xi_prim) * Math.sinh(2 * eta_prim) + beta2 * Math.cos(4 * xi_prim) * Math.sinh(4 * eta_prim) + beta3 * Math.cos(6 * xi_prim) * Math.sinh(6 * eta_prim) + beta4 * Math.cos(8 * xi_prim) * Math.sinh(8 * eta_prim)); return new EastNorth(y, x); } public LatLon eastNorth2latlon(EastNorth p) { double xi = (p.north() - false_n) / (k_zero * a_roof); double eta = (p.east() - false_e) / (k_zero * a_roof); double xi_prim = xi - delta1 * Math.sin(2 * xi) * Math.cosh(2 * eta) - delta2 * Math.sin(4 * xi) * Math.cosh(4 * eta) - delta3 * Math.sin(6 * xi) * Math.cosh(6 * eta) - delta4 * Math.sin(8 * xi) * Math.cosh(8 * eta); double eta_prim = eta - delta1 * Math.cos(2 * xi) * Math.sinh(2 * eta) - delta2 * Math.cos(4 * xi) * Math.sinh(4 * eta) - delta3 * Math.cos(6 * xi) * Math.sinh(6 * eta) - delta4 * Math.cos(8 * xi) * Math.sinh(8 * eta); double phi_star = Math.asin(Math.sin(xi_prim) / Math.cosh(eta_prim)); double delta_lambda = Math.atan(Math.sinh(eta_prim) / Math.cos(xi_prim)); double lambda = lambda_zero + delta_lambda; double phi = phi_star + Math.sin(phi_star) * Math.cos(phi_star) * (astar + bstar * Math.pow(Math.sin(phi_star), 2) + cstar * Math.pow(Math.sin(phi_star), 4) + dstar * Math.pow(Math.sin(phi_star), 6)); return new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda)); } public abstract String getCacheDirectoryName(); public abstract double scaleFactor(); private double atanh(double z) { return (Math.log(1 + z) - Math.log(1 - z)) / 2; } } package org.openstreetmap.josm.data.projection; public class Epsg3006 extends GaussKruger { public Epsg3006() { super(Ellipsoid.GRS80, 15.0, 0.0, 500000.0, 0.9996); } @Override public String getCacheDirectoryName() { return "Epsg:3006 (SWEREF 99 TM)"; } @Override public double scaleFactor() { /* Not sure what this method is supposed to do... The unit of Eastings is meters though. */ throw new UnsupportedOperationException("Not supported yet."); } } package org.openstreetmap.josm.data.projection; public class Epsg3007 extends GaussKruger { public Epsg3007() { super(Ellipsoid.GRS80, 12.0, 0.0, 500000.0, 1.0); } @Override public String getCacheDirectoryName() { return "Epsg:3007 (SWEREF 99 12 00)"; } @Override public double scaleFactor() { /* Not sure what this method is supposed to do... The unit of Eastings is meters though. */ throw new UnsupportedOperationException("Not supported yet."); } } _______________________________________________ Talk-se mailing list [email protected] http://lists.openstreetmap.org/listinfo/talk-se
