[4285] | 1 | // License: GPL. For details, see LICENSE file.
|
---|
| 2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
| 3 |
|
---|
[6362] | 4 | import static java.lang.Math.PI;
|
---|
| 5 | import static java.lang.Math.abs;
|
---|
| 6 | import static java.lang.Math.atan;
|
---|
| 7 | import static java.lang.Math.cos;
|
---|
| 8 | import static java.lang.Math.exp;
|
---|
| 9 | import static java.lang.Math.log;
|
---|
| 10 | import static java.lang.Math.pow;
|
---|
| 11 | import static java.lang.Math.sin;
|
---|
| 12 | import static java.lang.Math.sqrt;
|
---|
| 13 | import static java.lang.Math.tan;
|
---|
| 14 | import static java.lang.Math.toRadians;
|
---|
[4285] | 15 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
| 16 |
|
---|
[9124] | 17 | import org.openstreetmap.josm.data.Bounds;
|
---|
[8346] | 18 | import org.openstreetmap.josm.data.projection.CustomProjection.Param;
|
---|
[4285] | 19 | import org.openstreetmap.josm.data.projection.Ellipsoid;
|
---|
[5066] | 20 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
[4285] | 21 |
|
---|
| 22 | /**
|
---|
| 23 | * Implementation of the Lambert Conformal Conic projection.
|
---|
| 24 | *
|
---|
| 25 | * @author Pieren
|
---|
| 26 | */
|
---|
[9116] | 27 | public class LambertConformalConic extends AbstractProj {
|
---|
[5066] | 28 |
|
---|
[4285] | 29 | protected Ellipsoid ellps;
|
---|
[5066] | 30 |
|
---|
[6883] | 31 | public abstract static class Parameters {
|
---|
[5040] | 32 | public final double latitudeOrigin;
|
---|
[8510] | 33 |
|
---|
[5040] | 34 | public Parameters(double latitudeOrigin) {
|
---|
| 35 | this.latitudeOrigin = latitudeOrigin;
|
---|
| 36 | }
|
---|
[6142] | 37 | }
|
---|
[5066] | 38 |
|
---|
[5040] | 39 | public static class Parameters1SP extends Parameters {
|
---|
| 40 | public Parameters1SP(double latitudeOrigin) {
|
---|
| 41 | super(latitudeOrigin);
|
---|
| 42 | }
|
---|
| 43 | }
|
---|
[4285] | 44 |
|
---|
[5040] | 45 | public static class Parameters2SP extends Parameters {
|
---|
| 46 | public final double standardParallel1;
|
---|
| 47 | public final double standardParallel2;
|
---|
[8510] | 48 |
|
---|
[5040] | 49 | public Parameters2SP(double latitudeOrigin, double standardParallel1, double standardParallel2) {
|
---|
| 50 | super(latitudeOrigin);
|
---|
| 51 | this.standardParallel1 = standardParallel1;
|
---|
| 52 | this.standardParallel2 = standardParallel2;
|
---|
| 53 | }
|
---|
| 54 | }
|
---|
| 55 |
|
---|
| 56 | private Parameters params;
|
---|
[5066] | 57 |
|
---|
[4285] | 58 | /**
|
---|
| 59 | * projection exponent
|
---|
| 60 | */
|
---|
| 61 | protected double n;
|
---|
| 62 | /**
|
---|
| 63 | * projection factor
|
---|
| 64 | */
|
---|
[8346] | 65 | protected double f;
|
---|
[4285] | 66 | /**
|
---|
[5066] | 67 | * radius of the parallel of latitude of the false origin (2SP) or at
|
---|
[4285] | 68 | * natural origin (1SP)
|
---|
| 69 | */
|
---|
[5066] | 70 | protected double r0;
|
---|
| 71 |
|
---|
[4285] | 72 | /**
|
---|
| 73 | * precision in iterative schema
|
---|
| 74 | */
|
---|
| 75 | protected static final double epsilon = 1e-12;
|
---|
| 76 |
|
---|
[5066] | 77 | @Override
|
---|
| 78 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
[9600] | 79 | super.initialize(params);
|
---|
[5066] | 80 | ellps = params.ellps;
|
---|
[8346] | 81 | if (params.lat0 == null)
|
---|
| 82 | throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", Param.lat_0.key));
|
---|
| 83 | if (params.lat1 != null && params.lat2 != null) {
|
---|
| 84 | initialize2SP(params.lat0, params.lat1, params.lat2);
|
---|
[5066] | 85 | } else {
|
---|
[8346] | 86 | initialize1SP(params.lat0);
|
---|
[5066] | 87 | }
|
---|
[4285] | 88 | }
|
---|
| 89 |
|
---|
| 90 | /**
|
---|
| 91 | * Initialize for LCC with 2 standard parallels.
|
---|
[5066] | 92 | *
|
---|
[4285] | 93 | * @param lat_0 latitude of false origin (in degrees)
|
---|
| 94 | * @param lat_1 latitude of first standard parallel (in degrees)
|
---|
| 95 | * @param lat_2 latitude of second standard parallel (in degrees)
|
---|
| 96 | */
|
---|
[5066] | 97 | private void initialize2SP(double lat_0, double lat_1, double lat_2) {
|
---|
[5040] | 98 | this.params = new Parameters2SP(lat_0, lat_1, lat_2);
|
---|
[5066] | 99 |
|
---|
[4285] | 100 | final double m1 = m(toRadians(lat_1));
|
---|
| 101 | final double m2 = m(toRadians(lat_2));
|
---|
[5066] | 102 |
|
---|
[4285] | 103 | final double t1 = t(toRadians(lat_1));
|
---|
| 104 | final double t2 = t(toRadians(lat_2));
|
---|
| 105 | final double tf = t(toRadians(lat_0));
|
---|
[5066] | 106 |
|
---|
[4285] | 107 | n = (log(m1) - log(m2)) / (log(t1) - log(t2));
|
---|
[8346] | 108 | f = m1 / (n * pow(t1, n));
|
---|
| 109 | r0 = f * pow(tf, n);
|
---|
[4285] | 110 | }
|
---|
[5066] | 111 |
|
---|
[4285] | 112 | /**
|
---|
| 113 | * Initialize for LCC with 1 standard parallel.
|
---|
[5066] | 114 | *
|
---|
[4285] | 115 | * @param lat_0 latitude of natural origin (in degrees)
|
---|
| 116 | */
|
---|
[5066] | 117 | private void initialize1SP(double lat_0) {
|
---|
[5040] | 118 | this.params = new Parameters1SP(lat_0);
|
---|
[4285] | 119 | final double lat_0_rad = toRadians(lat_0);
|
---|
[5066] | 120 |
|
---|
[4285] | 121 | final double m0 = m(lat_0_rad);
|
---|
| 122 | final double t0 = t(lat_0_rad);
|
---|
[5066] | 123 |
|
---|
[4285] | 124 | n = sin(lat_0_rad);
|
---|
[8346] | 125 | f = m0 / (n * pow(t0, n));
|
---|
| 126 | r0 = f * pow(t0, n);
|
---|
[4285] | 127 | }
|
---|
| 128 |
|
---|
| 129 | /**
|
---|
| 130 | * auxiliary function t
|
---|
[9243] | 131 | * @param lat_rad latitude in radians
|
---|
| 132 | * @return result
|
---|
[4285] | 133 | */
|
---|
| 134 | protected double t(double lat_rad) {
|
---|
| 135 | return tan(PI/4 - lat_rad / 2.0)
|
---|
[8345] | 136 | / pow((1.0 - e * sin(lat_rad)) / (1.0 + e * sin(lat_rad)), e/2);
|
---|
[4285] | 137 | }
|
---|
| 138 |
|
---|
| 139 | /**
|
---|
| 140 | * auxiliary function m
|
---|
[9243] | 141 | * @param lat_rad latitude in radians
|
---|
| 142 | * @return result
|
---|
[4285] | 143 | */
|
---|
| 144 | protected double m(double lat_rad) {
|
---|
| 145 | return cos(lat_rad) / (sqrt(1 - e * e * pow(sin(lat_rad), 2)));
|
---|
| 146 | }
|
---|
[5066] | 147 |
|
---|
[4285] | 148 | @Override
|
---|
| 149 | public String getName() {
|
---|
| 150 | return tr("Lambert Conformal Conic");
|
---|
| 151 | }
|
---|
| 152 |
|
---|
| 153 | @Override
|
---|
| 154 | public String getProj4Id() {
|
---|
| 155 | return "lcc";
|
---|
| 156 | }
|
---|
| 157 |
|
---|
| 158 | @Override
|
---|
| 159 | public double[] project(double phi, double lambda) {
|
---|
| 160 | double sinphi = sin(phi);
|
---|
[8444] | 161 | double l = (0.5*log((1+sinphi)/(1-sinphi))) - e/2*log((1+e*sinphi)/(1-e*sinphi));
|
---|
| 162 | double r = f*exp(-n*l);
|
---|
[4285] | 163 | double gamma = n*lambda;
|
---|
[8444] | 164 | double x = r*sin(gamma);
|
---|
| 165 | double y = r0 - r*cos(gamma);
|
---|
| 166 | return new double[] {x, y};
|
---|
[4285] | 167 | }
|
---|
[5066] | 168 |
|
---|
[4285] | 169 | @Override
|
---|
| 170 | public double[] invproject(double east, double north) {
|
---|
[8510] | 171 | double r = sqrt(pow(east, 2) + pow(north-r0, 2));
|
---|
[4285] | 172 | double gamma = atan(east / (r0-north));
|
---|
| 173 | double lambda = gamma/n;
|
---|
[8346] | 174 | double latIso = (-1/n) * log(abs(r/f));
|
---|
[4285] | 175 | double phi = ellps.latitude(latIso, e, epsilon);
|
---|
[8444] | 176 | return new double[] {phi, lambda};
|
---|
[4285] | 177 | }
|
---|
[5066] | 178 |
|
---|
[5040] | 179 | public final Parameters getParameters() {
|
---|
| 180 | return params;
|
---|
| 181 | }
|
---|
[9124] | 182 |
|
---|
| 183 | @Override
|
---|
| 184 | public Bounds getAlgorithmBounds() {
|
---|
| 185 | double lat;
|
---|
| 186 | if (params instanceof Parameters2SP) {
|
---|
| 187 | Parameters2SP p2p = (Parameters2SP) params;
|
---|
| 188 | lat = (p2p.standardParallel1 + p2p.standardParallel2) / 2;
|
---|
| 189 | } else {
|
---|
| 190 | lat = params.latitudeOrigin;
|
---|
| 191 | }
|
---|
| 192 | double minlat = Math.max(lat - 60, -89);
|
---|
| 193 | double maxlat = Math.min(lat + 60, 89);
|
---|
| 194 | return new Bounds(minlat, -85, maxlat, 85, false);
|
---|
| 195 | }
|
---|
[4285] | 196 | }
|
---|