001// License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.data.projection.proj; 003 004import static java.lang.Math.PI; 005import static java.lang.Math.abs; 006import static java.lang.Math.asin; 007import static java.lang.Math.atan; 008import static java.lang.Math.atan2; 009import static java.lang.Math.cos; 010import static java.lang.Math.exp; 011import static java.lang.Math.log; 012import static java.lang.Math.pow; 013import static java.lang.Math.sin; 014import static java.lang.Math.sqrt; 015import static java.lang.Math.tan; 016import static org.openstreetmap.josm.tools.I18n.tr; 017 018import org.openstreetmap.josm.data.Bounds; 019import org.openstreetmap.josm.data.projection.Ellipsoid; 020import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; 021import org.openstreetmap.josm.tools.JosmRuntimeException; 022import org.openstreetmap.josm.tools.Utils; 023 024// CHECKSTYLE.OFF: LineLength 025 026/** 027 * Projection for the SwissGrid CH1903 / L03, see <a href="https://en.wikipedia.org/wiki/Swiss_coordinate_system">Wikipedia article</a>.<br> 028 * 029 * Calculations were originally based on 030 * <a href="http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf"> 031 * simple formula</a>.<br> 032 * 033 * August 2010 update to 034 * <a href="http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf"> 035 * this formula (rigorous formulas)</a>. 036 */ 037public class SwissObliqueMercator extends AbstractProj { 038 039 // CHECKSTYLE.ON: LineLength 040 041 private Ellipsoid ellps; 042 private double kR; 043 private double alpha; 044 private double b0; 045 private double k; 046 private double phi0; 047 048 private static final double EPSILON = 1e-11; 049 050 @Override 051 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 052 super.initialize(params); 053 if (params.lat0 == null) 054 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 055 ellps = params.ellps; 056 initialize(params.lat0); 057 } 058 059 private void initialize(double lat0) { 060 phi0 = Utils.toRadians(lat0); 061 kR = sqrt(1 - ellps.e2) / (1 - (ellps.e2 * pow(sin(phi0), 2))); 062 alpha = sqrt(1 + (ellps.eb2 * pow(cos(phi0), 4))); 063 b0 = asin(sin(phi0) / alpha); 064 k = log(tan(PI / 4 + b0 / 2)) - alpha 065 * log(tan(PI / 4 + phi0 / 2)) + alpha * ellps.e / 2 066 * log((1 + ellps.e * sin(phi0)) / (1 - ellps.e * sin(phi0))); 067 } 068 069 @Override 070 public String getName() { 071 return tr("Swiss Oblique Mercator"); 072 } 073 074 @Override 075 public String getProj4Id() { 076 return "somerc"; 077 } 078 079 @Override 080 public double[] project(double phi, double lambda) { 081 double s = alpha * log(tan(PI / 4 + phi / 2)) - alpha * ellps.e / 2 082 * log((1 + ellps.e * sin(phi)) / (1 - ellps.e * sin(phi))) + k; 083 double b = 2 * (atan(exp(s)) - PI / 4); 084 double l = alpha * lambda; 085 086 double lb = atan2(sin(l), sin(b0) * tan(b) + cos(b0) * cos(l)); 087 double bb = asin(cos(b0) * sin(b) - sin(b0) * cos(b) * cos(l)); 088 089 double y = kR * lb; 090 double x = kR / 2 * log((1 + sin(bb)) / (1 - sin(bb))); 091 092 return new double[] {y, x}; 093 } 094 095 @Override 096 public double[] invproject(double y, double x) { 097 double lb = y / kR; 098 double bb = 2 * (atan(exp(x / kR)) - PI / 4); 099 100 double b = asin(cos(b0) * sin(bb) + sin(b0) * cos(bb) * cos(lb)); 101 double l = atan2(sin(lb), cos(b0) * cos(lb) - sin(b0) * tan(bb)); 102 103 double lambda = l / alpha; 104 double phi = b; 105 106 double prevPhi = -1000; 107 int iteration = 0; 108 // iteration to finds S and phi 109 while (abs(phi - prevPhi) > EPSILON) { 110 if (++iteration > 30) 111 throw new JosmRuntimeException("Too many iterations"); 112 prevPhi = phi; 113 double s = 1 / alpha * (log(tan(PI / 4 + b / 2)) - k) + ellps.e 114 * log(tan(PI / 4 + asin(ellps.e * sin(phi)) / 2)); 115 phi = 2 * atan(exp(s)) - PI / 2; 116 } 117 return new double[] {phi, lambda}; 118 } 119 120 @Override 121 public Bounds getAlgorithmBounds() { 122 if (phi0 > 0) { 123 return new Bounds(-10, -40, 85, 40, false); 124 } else { 125 return new Bounds(-85, -40, 10, 40, false); 126 } 127 } 128}