| 1 | //License: GPL. For details, see LICENSE file.
|
|---|
| 2 | //Thanks to Johan Montagnat and its geoconv java converter application
|
|---|
| 3 | //(http://www.i3s.unice.fr/~johan/gps/ , published under GPL license)
|
|---|
| 4 | //from which some code and constants have been reused here.
|
|---|
| 5 | package org.openstreetmap.josm.data.projection;
|
|---|
| 6 |
|
|---|
| 7 | import static org.openstreetmap.josm.tools.I18n.tr;
|
|---|
| 8 |
|
|---|
| 9 | import org.openstreetmap.josm.data.coor.EastNorth;
|
|---|
| 10 | import org.openstreetmap.josm.data.coor.LatLon;
|
|---|
| 11 |
|
|---|
| 12 | public class LambertEST implements Projection {
|
|---|
| 13 |
|
|---|
| 14 | public static final double ef = 500000; //Easting at false origin = 5000000 m
|
|---|
| 15 | public static final double nf = 6375000; //Northing at false origin = 6375000 m
|
|---|
| 16 | public static final double lat1 = Math.toRadians(59 + 1.0/3.0); //Latitude of 1st standard parallel = 59o20`0`` N
|
|---|
| 17 | public static final double lat2 = Math.toRadians( 58);//Latitude of 2nd standard parallel = 58o0`0`` N
|
|---|
| 18 | public static final double latf = Math.toRadians(57.517553930555555555555555555556);//'Latitude of false origin = 57o31`3,19415`` N
|
|---|
| 19 | public static final double lonf = Math.toRadians( 24.0);
|
|---|
| 20 | public static final double a = 6378137;
|
|---|
| 21 | public static final double ee = 0.081819191042815792;
|
|---|
| 22 | public static final double m1 = Math.cos(lat1) / (Math.sqrt(1 - ee *ee * Math.pow(Math.sin(lat1), 2)));
|
|---|
| 23 | public static final double m2 = Math.cos(lat2) / (Math.sqrt(1 - ee *ee * Math.pow(Math.sin(lat2), 2)));
|
|---|
| 24 | public static final double t1 = Math.tan(Math.PI / 4.0 - lat1 / 2.0) / Math.pow(( (1.0 - ee * Math.sin(lat1)) / (1.0 + ee * Math.sin(lat1))) ,(ee / 2.0));
|
|---|
| 25 | public static final double t2 = Math.tan(Math.PI / 4.0 - lat2 / 2.0) / Math.pow(( (1.0 - ee * Math.sin(lat2)) / (1.0 + ee * Math.sin(lat2))) ,(ee / 2.0));
|
|---|
| 26 | public static final double tf = Math.tan(Math.PI / 4.0 - latf / 2.0) / Math.pow(( (1.0 - ee * Math.sin(latf)) / (1.0 + ee * Math.sin(latf))) ,(ee / 2.0));
|
|---|
| 27 | public static final double n = (Math.log(m1) - Math.log(m2)) / (Math.log(t1) - Math.log(t2));
|
|---|
| 28 | public static final double f = m1 / (n * Math.pow(t1, n));
|
|---|
| 29 | public static final double rf = a * f * Math.pow(tf, n);
|
|---|
| 30 |
|
|---|
| 31 | // /**
|
|---|
| 32 | // * precision in iterative schema
|
|---|
| 33 | // */
|
|---|
| 34 | //
|
|---|
| 35 | public static final double epsilon = 1e-11;
|
|---|
| 36 |
|
|---|
| 37 | /**
|
|---|
| 38 | * @param p WGS84 lat/lon (ellipsoid GRS80) (in degree)
|
|---|
| 39 | * @return eastnorth projection in Lambert Zone (ellipsoid GRS80)
|
|---|
| 40 | */
|
|---|
| 41 | public EastNorth latlon2eastNorth(LatLon p)
|
|---|
| 42 | {
|
|---|
| 43 |
|
|---|
| 44 | double t = Math.tan(Math.PI / 4.0 - Math.toRadians(p.lat()) / 2.0) / Math.pow(( (1.0 - ee * Math.sin(Math.toRadians(p.lat()))) / (1.0 + ee * Math.sin(Math.toRadians(p.lat())))) ,(ee / 2.0));
|
|---|
| 45 | double r = a * f * Math.pow(t, n);
|
|---|
| 46 | double theta = n * (Math.toRadians(p.lon()) - lonf);
|
|---|
| 47 |
|
|---|
| 48 | double x = ef + r * Math.sin(theta); //587446.7
|
|---|
| 49 | double y = nf + rf - r * Math.cos(theta); //6485401.6
|
|---|
| 50 |
|
|---|
| 51 | return new EastNorth(x,y);
|
|---|
| 52 | }
|
|---|
| 53 |
|
|---|
| 54 | public static double IterateAngle(double e, double t)
|
|---|
| 55 | {
|
|---|
| 56 | double a1 = 0.0;
|
|---|
| 57 | double a2 = 3.1415926535897931;
|
|---|
| 58 | double a = 1.5707963267948966;
|
|---|
| 59 | double b = 1.5707963267948966 - (2.0 * Math.atan(t * Math.pow((1.0 - (e * Math.sin(a))) / (1.0 + (e * Math.sin(a))), e / 2.0)));
|
|---|
| 60 | while (Math.abs(a-b) > epsilon)
|
|---|
| 61 | {
|
|---|
| 62 | a = a1 + ((a2 - a1) / 2.0);
|
|---|
| 63 | b = 1.5707963267948966 - (2.0 * Math.atan(t * Math.pow((1.0 - (e * Math.sin(a))) / (1.0 + (e * Math.sin(a))), e / 2.0)));
|
|---|
| 64 | if (a1 == a2)
|
|---|
| 65 | {
|
|---|
| 66 | return 0.0;
|
|---|
| 67 | }
|
|---|
| 68 | if (b > a)
|
|---|
| 69 | a1 = a;
|
|---|
| 70 | else
|
|---|
| 71 | a2 = a;
|
|---|
| 72 | }
|
|---|
| 73 | return b;
|
|---|
| 74 | }
|
|---|
| 75 | public LatLon eastNorth2latlon(EastNorth p)
|
|---|
| 76 | {
|
|---|
| 77 | double r = Math.sqrt(Math.pow((p.getX() - ef), 2.0) + Math.pow((rf - p.getY() + nf), 2.0) ) * Math.signum(n);
|
|---|
| 78 | double T = Math.pow((r / (a * f)), (1.0/ n)) ;
|
|---|
| 79 | double theta = Math.atan((p.getX() - ef) / (rf - p.getY() + nf));
|
|---|
| 80 | double y = (theta / n + lonf) ;
|
|---|
| 81 | double x = (IterateAngle(ee, T)) ;
|
|---|
| 82 | return new LatLon(Math.toDegrees(x),Math.toDegrees(y));
|
|---|
| 83 | }
|
|---|
| 84 |
|
|---|
| 85 | @Override public String toString() {
|
|---|
| 86 | return tr("Lambert Zone (Estonia)");
|
|---|
| 87 | }
|
|---|
| 88 |
|
|---|
| 89 | public String getCacheDirectoryName() {
|
|---|
| 90 | return "lambertest";
|
|---|
| 91 | }
|
|---|
| 92 |
|
|---|
| 93 | public double scaleFactor() {
|
|---|
| 94 | return 1.0 / 360;
|
|---|
| 95 | }
|
|---|
| 96 |
|
|---|
| 97 | @Override
|
|---|
| 98 | public boolean equals(Object o) {
|
|---|
| 99 | return o instanceof LambertEST;
|
|---|
| 100 | }
|
|---|
| 101 |
|
|---|
| 102 | @Override
|
|---|
| 103 | public int hashCode() {
|
|---|
| 104 | return LambertEST.class.hashCode();
|
|---|
| 105 | }
|
|---|
| 106 | }
|
|---|