001 // License: GPL. For details, see LICENSE file. 002 package org.openstreetmap.josm.data.projection.proj; 003 004 import static java.lang.Math.*; 005 006 import static org.openstreetmap.josm.tools.I18n.tr; 007 008 import org.openstreetmap.josm.data.projection.Ellipsoid; 009 import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; 010 011 /** 012 * Implementation of the Lambert Conformal Conic projection. 013 * 014 * @author Pieren 015 */ 016 public class LambertConformalConic implements Proj { 017 018 protected Ellipsoid ellps; 019 protected double e; 020 021 public static abstract class Parameters { 022 public final double latitudeOrigin; 023 public Parameters(double latitudeOrigin) { 024 this.latitudeOrigin = latitudeOrigin; 025 } 026 }; 027 028 public static class Parameters1SP extends Parameters { 029 public Parameters1SP(double latitudeOrigin) { 030 super(latitudeOrigin); 031 } 032 } 033 034 public static class Parameters2SP extends Parameters { 035 public final double standardParallel1; 036 public final double standardParallel2; 037 public Parameters2SP(double latitudeOrigin, double standardParallel1, double standardParallel2) { 038 super(latitudeOrigin); 039 this.standardParallel1 = standardParallel1; 040 this.standardParallel2 = standardParallel2; 041 } 042 } 043 044 private Parameters params; 045 046 /** 047 * projection exponent 048 */ 049 protected double n; 050 /** 051 * projection factor 052 */ 053 protected double F; 054 /** 055 * radius of the parallel of latitude of the false origin (2SP) or at 056 * natural origin (1SP) 057 */ 058 protected double r0; 059 060 /** 061 * precision in iterative schema 062 */ 063 protected static final double epsilon = 1e-12; 064 065 @Override 066 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 067 ellps = params.ellps; 068 e = ellps.e; 069 if (params.lat_0 == null) 070 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 071 if (params.lat_1 != null && params.lat_2 != null) { 072 initialize2SP(params.lat_0, params.lat_1, params.lat_2); 073 } else { 074 initialize1SP(params.lat_0); 075 } 076 } 077 078 /** 079 * Initialize for LCC with 2 standard parallels. 080 * 081 * @param lat_0 latitude of false origin (in degrees) 082 * @param lat_1 latitude of first standard parallel (in degrees) 083 * @param lat_2 latitude of second standard parallel (in degrees) 084 */ 085 private void initialize2SP(double lat_0, double lat_1, double lat_2) { 086 this.params = new Parameters2SP(lat_0, lat_1, lat_2); 087 088 final double m1 = m(toRadians(lat_1)); 089 final double m2 = m(toRadians(lat_2)); 090 091 final double t1 = t(toRadians(lat_1)); 092 final double t2 = t(toRadians(lat_2)); 093 final double tf = t(toRadians(lat_0)); 094 095 n = (log(m1) - log(m2)) / (log(t1) - log(t2)); 096 F = m1 / (n * pow(t1, n)); 097 r0 = F * pow(tf, n); 098 } 099 100 /** 101 * Initialize for LCC with 1 standard parallel. 102 * 103 * @param lat_0 latitude of natural origin (in degrees) 104 */ 105 private void initialize1SP(double lat_0) { 106 this.params = new Parameters1SP(lat_0); 107 final double lat_0_rad = toRadians(lat_0); 108 109 final double m0 = m(lat_0_rad); 110 final double t0 = t(lat_0_rad); 111 112 n = sin(lat_0_rad); 113 F = m0 / (n * pow(t0, n)); 114 r0 = F * pow(t0, n); 115 } 116 117 /** 118 * auxiliary function t 119 */ 120 protected double t(double lat_rad) { 121 return tan(PI/4 - lat_rad / 2.0) 122 / pow(( (1.0 - e * sin(lat_rad)) / (1.0 + e * sin(lat_rad))) , e/2); 123 } 124 125 /** 126 * auxiliary function m 127 */ 128 protected double m(double lat_rad) { 129 return cos(lat_rad) / (sqrt(1 - e * e * pow(sin(lat_rad), 2))); 130 } 131 132 @Override 133 public String getName() { 134 return tr("Lambert Conformal Conic"); 135 } 136 137 @Override 138 public String getProj4Id() { 139 return "lcc"; 140 } 141 142 @Override 143 public double[] project(double phi, double lambda) { 144 double sinphi = sin(phi); 145 double L = (0.5*log((1+sinphi)/(1-sinphi))) - e/2*log((1+e*sinphi)/(1-e*sinphi)); 146 double r = F*exp(-n*L); 147 double gamma = n*lambda; 148 double X = r*sin(gamma); 149 double Y = r0 - r*cos(gamma); 150 return new double[] { X, Y }; 151 } 152 153 @Override 154 public double[] invproject(double east, double north) { 155 double r = sqrt(pow(east,2) + pow(north-r0, 2)); 156 double gamma = atan(east / (r0-north)); 157 double lambda = gamma/n; 158 double latIso = (-1/n) * log(abs(r/F)); 159 double phi = ellps.latitude(latIso, e, epsilon); 160 return new double[] { phi, lambda }; 161 } 162 163 public final Parameters getParameters() { 164 return params; 165 } 166 }