Index: trunk/src/org/openstreetmap/josm/data/projection/AbstractProjection.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/AbstractProjection.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/AbstractProjection.java (revision 4285)
@@ -0,0 +1,83 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection;
+
+import org.openstreetmap.josm.data.projection.datum.Datum;
+import org.openstreetmap.josm.data.projection.proj.Proj;
+import org.openstreetmap.josm.data.coor.EastNorth;
+import org.openstreetmap.josm.data.coor.LatLon;
+
+/**
+ * Implementation of the Projection interface that represents a coordinate reference system and delegates
+ * the real projection and datum conversion to other classes.
+ *
+ * It handles false easting and northing, central meridian and general scale factor before calling the
+ * delegate projection.
+ *
+ * Forwards lat/lon values to the real projection in units of radians.
+ *
+ * The fields are named after Proj.4 parameters.
+ *
+ * Subclasses of AbstractProjection must set ellps and proj to a non-null value.
+ * In addition, either datum or nadgrid has to be initialized to some value.
+ *
+ * FIXME: nadgrids should probably be implemented as a Datum
+ */
+abstract public class AbstractProjection implements Projection {
+
+ protected Ellipsoid ellps;
+ protected Datum datum;
+ protected Proj proj;
+ protected double x_0 = 0.0; /* false easting (in meters) */
+ protected double y_0 = 0.0; /* false northing (in meters) */
+ protected double lon_0 = 0.0; /* central meridian */
+ protected double k_0 = 1.0; /* general scale factor */
+ protected NTV2GridShiftFile nadgrids = null;
+
+ @Override
+ public EastNorth latlon2eastNorth(LatLon ll) {
+ if (nadgrids != null) {
+ NTV2GridShift gs = new NTV2GridShift(ll);
+ nadgrids.gridShiftReverse(gs);
+ ll = new LatLon(ll.lat()+gs.getLatShiftDegrees(), ll.lon()+gs.getLonShiftPositiveEastDegrees());
+ } else {
+ ll = datum.fromWGS84(ll);
+ }
+ double[] en = proj.project(Math.toRadians(ll.lat()), Math.toRadians(ll.lon() - lon_0));
+ return new EastNorth(ellps.a * k_0 * en[0] + x_0, ellps.a * k_0 * en[1] + y_0);
+ }
+
+ @Override
+ public LatLon eastNorth2latlon(EastNorth en) {
+ double[] latlon_rad = proj.invproject((en.east() - x_0) / ellps.a / k_0, (en.north() - y_0) / ellps.a / k_0);
+ LatLon ll = new LatLon(Math.toDegrees(latlon_rad[0]), Math.toDegrees(latlon_rad[1]) + lon_0);
+ if (nadgrids != null) {
+ NTV2GridShift gs = new NTV2GridShift(ll);
+ nadgrids.gridShiftForward(gs);
+ ll = new LatLon(ll.lat()+gs.getLatShiftDegrees(), ll.lon()+gs.getLonShiftPositiveEastDegrees());
+ } else {
+ ll = datum.toWGS84(ll);
+ }
+ return ll;
+ }
+
+ @Override
+ public double getDefaultZoomInPPD() {
+ // this will set the map scaler to about 1000 m
+ return 10;
+ }
+
+ /**
+ * @return The EPSG Code of this CRS, null if it doesn't have one.
+ */
+ public abstract Integer getEpsgCode();
+
+ /**
+ * Default implementation of toCode().
+ * Should be overridden, if there is no EPSG code for this CRS.
+ */
+ @Override
+ public String toCode() {
+ return "EPSG:" + getEpsgCode();
+ }
+
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java (revision 4285)
@@ -19,5 +19,6 @@
public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
/**
- * Hayford's ellipsoid (ED50 system)
+ * Hayford's ellipsoid 1909 (ED50 system)
+ * Proj.4 code: intl
*/
public static final Ellipsoid hayford =
Index: trunk/src/org/openstreetmap/josm/data/projection/Epsg3008.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/Epsg3008.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/Epsg3008.java (revision 4285)
@@ -6,5 +6,6 @@
import org.openstreetmap.josm.data.Bounds;
import org.openstreetmap.josm.data.coor.LatLon;
-import org.openstreetmap.josm.data.coor.EastNorth;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
+import org.openstreetmap.josm.data.projection.proj.TransverseMercator;
/**
@@ -12,29 +13,24 @@
* http://spatialreference.org/ref/epsg/3008/
*
- * @author Hanno Hecker, based on the TransverseMercatorLV.java by Viesturs Zarins
+ * @author Hanno Hecker
*/
-public class Epsg3008 extends TransverseMercator {
+public class Epsg3008 extends AbstractProjection {
- private final static double UTMScaleFactor = 1.0;
- private double UTMCentralMeridianRad;
- private double offsetEastMeters = 150000;
- private double offsetNorthMeters = 0;
-
- public Epsg3008()
- {
- UTMCentralMeridianRad = Math.toRadians(13.5);
+ public Epsg3008() {
+ ellps = Ellipsoid.GRS80;
+ proj = new TransverseMercator(ellps);
+ datum = GRS80Datum.INSTANCE;
+ lon_0 = 13.5;
+ x_0 = 150000;
}
-
- @Override public String toString() {
+
+ @Override
+ public String toString() {
return tr("SWEREF99 13 30 / EPSG:3008 (Sweden)");
}
- private int epsgCode() {
+ @Override
+ public Integer getEpsgCode() {
return 3008;
- }
-
- @Override
- public String toCode() {
- return "EPSG:"+ epsgCode();
}
@@ -44,6 +40,7 @@
}
+ @Override
public String getCacheDirectoryName() {
- return "epsg"+ epsgCode();
+ return "epsg"+ getEpsgCode();
}
@@ -55,14 +52,3 @@
}
- @Override
- public EastNorth latlon2eastNorth(LatLon p) {
- EastNorth a = mapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), UTMCentralMeridianRad);
- return new EastNorth(a.east() * UTMScaleFactor + offsetEastMeters, a.north() * UTMScaleFactor + offsetNorthMeters);
- }
-
- @Override
- public LatLon eastNorth2latlon(EastNorth p) {
- return mapXYToLatLon((p.east() - offsetEastMeters)/UTMScaleFactor, (p.north() - offsetNorthMeters)/UTMScaleFactor, UTMCentralMeridianRad);
- }
-
}
Index: trunk/src/org/openstreetmap/josm/data/projection/Lambert.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/Lambert.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/Lambert.java (revision 4285)
@@ -6,5 +6,4 @@
import java.awt.GridBagLayout;
import java.awt.event.ActionListener;
-import java.io.IOException;
import java.io.InputStream;
import java.util.Collection;
@@ -17,48 +16,35 @@
import org.openstreetmap.josm.Main;
import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.proj.LambertConformalConic;
import org.openstreetmap.josm.tools.GBC;
import org.openstreetmap.josm.tools.ImageProvider;
/**
- * This class provides the two methods latlon2eastNorth and eastNorth2latlon
- * converting the JOSM LatLon coordinates in WGS84 system (GPS) to and from East North values in
- * the projection Lambert conic conform 4 zones using the French geodetic system NTF.
+ * Lambert conic conform 4 zones using the French geodetic system NTF.
* This newer version uses the grid translation NTF<->RGF93 provided by IGN for a submillimetric accuracy.
* (RGF93 is the French geodetic system similar to WGS84 but not mathematically equal)
* @author Pieren
*/
-public class Lambert implements Projection, ProjectionSubPrefs {
+public class Lambert extends AbstractProjection implements ProjectionSubPrefs {
/**
* Lambert I, II, III, and IV projection exponents
*/
- public static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 };
+ private static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 };
/**
* Lambert I, II, III, and IV projection constants
*/
- public static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 };
+ private static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 };
/**
* Lambert I, II, III, and IV false east
*/
- public static final double Xs[] = { 600000.0, 600000.0, 600000.0, 234.358 };
+ private static final double x_0s[] = { 600000.0, 600000.0, 600000.0, 234.358 };
/**
* Lambert I, II, III, and IV false north
*/
- public static final double Ys[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 };
-
- /**
- * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian
- */
- public static final double lg0 = 0.04079234433198; // 2deg20'14.025"
-
- /**
- * precision in iterative schema
- */
-
- public static final double epsilon = 1e-11;
+ private static final double y_fs[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 };
/**
@@ -68,5 +54,5 @@
public static final double cMinLatZone1Radian = Math.toRadians(46.1 * 0.9);// lowest latitude of Zone 4 (South Corsica)
- public static final double zoneLimitsDegree[][] = {
+ public static final double[][] zoneLimitsDegree = {
{Math.toDegrees(cMaxLatZone1Radian), (53.5 * 0.9)}, // Zone 1 (reference values in grad *0.9)
{(53.5 * 0.9), (50.5 * 0.9)}, // Zone 2
@@ -86,5 +72,5 @@
public static final int DEFAULT_ZONE = 0;
- private int layoutZone = DEFAULT_ZONE;
+ private int layoutZone;
private static NTV2GridShiftFile ntf_rgf93Grid = null;
@@ -100,46 +86,37 @@
InputStream is = Main.class.getResourceAsStream("/data/"+gridFileName);
if (is == null) {
- System.err.println(tr("Warning: failed to open input stream for resource ''/data/{0}''. Cannot load NTF<->RGF93 grid", gridFileName));
- return;
+ throw new RuntimeException(tr("Warning: failed to open input stream for resource ''/data/{0}''. Cannot load NTF<->RGF93 grid", gridFileName));
}
ntf_rgf93Grid = new NTV2GridShiftFile();
ntf_rgf93Grid.loadGridShiftFile(is, false);
- //System.out.println("NTF<->RGF93 grid loaded.");
} catch (Exception e) {
- e.printStackTrace();
+ throw new RuntimeException(e);
}
}
- }
-
- /**
- * @param p WGS84 lat/lon (ellipsoid GRS80) (in degree)
- * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
- * @throws IOException
- */
- public EastNorth latlon2eastNorth(LatLon p) {
- // translate ellipsoid GRS80 (WGS83) => Clark
- LatLon geo = WGS84_to_NTF(p);
- double lt = Math.toRadians(geo.lat()); // in radian
- double lg = Math.toRadians(geo.lon());
-
- // check if longitude and latitude are inside the French Lambert zones
- if (lt >= cMinLatZone1Radian && lt <= cMaxLatZone1Radian && lg >= cMinLonZonesRadian && lg <= cMaxLonZonesRadian)
- return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
- return ConicProjection(lt, lg, Xs[0], Ys[0], c[0], n[0]);
- }
-
- public LatLon eastNorth2latlon(EastNorth p) {
- LatLon geo;
- geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
- // translate geodetic system from NTF (ellipsoid Clark) to RGF93/WGS84 (ellipsoid GRS80)
- return NTF_to_WGS84(geo);
- }
-
- @Override public String toString() {
+ updateParameters(DEFAULT_ZONE);
+ }
+
+ private void updateParameters(int layoutZone) {
+ this.layoutZone = layoutZone;
+ ellps = Ellipsoid.clarke;
+ datum = null; // no datum needed, we have a shift file
+ nadgrids = ntf_rgf93Grid;
+ x_0 = x_0s[layoutZone];
+ lon_0 = 2.0 + 20.0 / 60 + 14.025 / 3600; // 0 grade Paris
+ if (proj == null) {
+ proj = new LambertConformalConic();
+ }
+ ((LambertConformalConic)proj).updateParametersDirect(
+ Ellipsoid.clarke, n[layoutZone], c[layoutZone] / ellps.a, y_fs[layoutZone] / ellps.a);
+ }
+
+ @Override
+ public String toString() {
return tr("Lambert 4 Zones (France)");
}
- public String toCode() {
- return "EPSG:"+(27561+layoutZone);
+ @Override
+ public Integer getEpsgCode() {
+ return 27561+layoutZone;
}
@@ -149,83 +126,10 @@
}
+ @Override
public String getCacheDirectoryName() {
return "lambert";
}
- /**
- * Initializes from geographic coordinates. Note that reference ellipsoid
- * used by Lambert is the Clark ellipsoid.
- *
- * @param lat latitude in grad
- * @param lon longitude in grad
- * @param Xs false east (coordinate system origin) in meters
- * @param Ys false north (coordinate system origin) in meters
- * @param c projection constant
- * @param n projection exponent
- * @return EastNorth projected coordinates in meter
- */
- private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) {
- double eslt = Ellipsoid.clarke.e * Math.sin(lat);
- double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0))
- * Math.pow((1.0 - eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0));
- double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0));
- double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0));
- return new EastNorth(east, north);
- }
-
- /**
- * Initializes from projected coordinates (conic projection). Note that
- * reference ellipsoid used by Lambert is Clark
- *
- * @param eastNorth projected coordinates pair in meters
- * @param Xs false east (coordinate system origin) in meters
- * @param Ys false north (coordinate system origin) in meters
- * @param c projection constant
- * @param n projection exponent
- * @return LatLon in degree
- */
- private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) {
- double dx = eastNorth.east() - Xs;
- double dy = Ys - eastNorth.north();
- double R = Math.sqrt(dx * dx + dy * dy);
- double gamma = Math.atan(dx / dy);
- double l = -1.0 / n * Math.log(Math.abs(R / c));
- l = Math.exp(l);
- double lon = lg0 + gamma / n;
- double lat = 2.0 * Math.atan(l) - Math.PI / 2.0;
- double delta = 1.0;
- while (delta > epsilon) {
- double eslt = Ellipsoid.clarke.e * Math.sin(lat);
- double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0 - eslt), Ellipsoid.clarke.e / 2.0) * l) - Math.PI
- / 2.0;
- delta = Math.abs(nlt - lat);
- lat = nlt;
- }
- return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon)); // in radian
- }
-
- /**
- * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert
- * geographic, (ellipsoid Clark)
- * @throws IOException
- */
- private LatLon WGS84_to_NTF(LatLon wgs) {
- NTV2GridShift gs = new NTV2GridShift(wgs);
- if (ntf_rgf93Grid != null) {
- ntf_rgf93Grid.gridShiftReverse(gs);
- return new LatLon(wgs.lat()+gs.getLatShiftDegrees(), wgs.lon()+gs.getLonShiftPositiveEastDegrees());
- }
- return new LatLon(0,0);
- }
-
- private LatLon NTF_to_WGS84(LatLon ntf) {
- NTV2GridShift gs = new NTV2GridShift(ntf);
- if (ntf_rgf93Grid != null) {
- ntf_rgf93Grid.gridShiftForward(gs);
- return new LatLon(ntf.lat()+gs.getLatShiftDegrees(), ntf.lon()+gs.getLonShiftPositiveEastDegrees());
- }
- return new LatLon(0,0);
- }
-
+ @Override
public Bounds getWorldBoundsLatLon()
{
@@ -234,12 +138,4 @@
new LatLon(Math.min(zoneLimitsDegree[layoutZone][0] + cMaxOverlappingZonesDegree, Math.toDegrees(cMaxLatZone1Radian)), Math.toDegrees(cMaxLonZonesRadian)));
return b;
- }
-
- /**
- * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
- */
- public double getDefaultZoomInPPD() {
- // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
- return 10.0;
}
@@ -273,4 +169,5 @@
}
+ @Override
public Collection getPreferences(JPanel p) {
Object prefcb = p.getComponent(2);
@@ -281,6 +178,7 @@
}
+ @Override
public void setPreferences(Collection args) {
- layoutZone = DEFAULT_ZONE;
+ int layoutZone = DEFAULT_ZONE;
if (args != null) {
try {
@@ -295,4 +193,5 @@
} catch(NumberFormatException e) {}
}
+ updateParameters(layoutZone);
}
@@ -306,4 +205,5 @@
}
+ @Override
public Collection getPreferencesFromCode(String code) {
if (code.startsWith("EPSG:2756") && code.length() == 10) {
Index: trunk/src/org/openstreetmap/josm/data/projection/LambertCC9Zones.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/LambertCC9Zones.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/LambertCC9Zones.java (revision 4285)
@@ -14,50 +14,17 @@
import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
+import org.openstreetmap.josm.data.projection.proj.LambertConformalConic;
import org.openstreetmap.josm.tools.GBC;
import org.openstreetmap.josm.tools.ImageProvider;
/**
- * This class implements the Lambert Conic Conform 9 Zones projection as specified by the IGN
+ * Lambert Conic Conform 9 Zones projection as specified by the IGN
* in this document http://professionnels.ign.fr/DISPLAY/000/526/700/5267002/transformation.pdf
* @author Pieren
*
*/
-public class LambertCC9Zones implements Projection, ProjectionSubPrefs {
-
- /**
- * Lambert 9 zones projection exponents
- */
- public static final double n[] = { 0.6691500006885269, 0.682018118346418, 0.6946784863203991, 0.7071272481559119,
- 0.7193606118567315, 0.7313748510399917, 0.7431663060711892, 0.7547313851789208, 0.7660665655489937};
-
- /**
- * Lambert 9 zones projection constants
- */
- public static final double c[] = { 1.215363305807804E7, 1.2050261119223533E7, 1.195716926884592E7, 1.18737533925172E7,
- 1.1799460698022118E7, 1.17337838820243E7, 1.16762559948139E7, 1.1626445901183508E7, 1.1583954251630554E7};
-
- /**
- * Lambert 9 zones false east
- */
- public static final double Xs = 1700000;
-
- /**
- * Lambert 9 zones false north
- */
- public static final double Ys[] = { 8293467.503439436, 9049604.665107645, 9814691.693461388, 1.0588107871787189E7,
- 1.1369285637569271E7, 1.2157704903382052E7, 1.2952888086405803E7, 1.3754395745267643E7, 1.4561822739114787E7};
-
- /**
- * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian
- */
- public static final double lg0 = 0.04079234433198; // 2deg20'14.025"
-
- /**
- * precision in iterative schema
- */
-
- public static final double epsilon = 1e-12;
+public class LambertCC9Zones extends AbstractProjection implements ProjectionSubPrefs {
/**
@@ -67,14 +34,4 @@
public static final double cMinLatZonesDegree = 41.0;
- public static final double cMinLatZonesRadian = Math.toRadians(cMinLatZonesDegree);
-
- public static final double cMinLonZonesRadian = Math.toRadians(-5.0);
-
- public static final double cMaxLonZonesRadian = Math.toRadians(10.2);
-
- public static final double lambda0 = Math.toRadians(3);
- public static final double e = Ellipsoid.GRS80.e; // but in doc=0.08181919112
- public static final double e2 =Ellipsoid.GRS80.e2;
- public static final double a = Ellipsoid.GRS80.a;
public static final double cMaxOverlappingZones = 1.5;
@@ -82,52 +39,32 @@
public static final int DEFAULT_ZONE = 0;
- private static int layoutZone = DEFAULT_ZONE;
+ private int layoutZone = DEFAULT_ZONE;
- private double L(double phi, double e) {
- double sinphi = Math.sin(phi);
- return (0.5*Math.log((1+sinphi)/(1-sinphi))) - e/2*Math.log((1+e*sinphi)/(1-e*sinphi));
+ public LambertCC9Zones() {
+ this(DEFAULT_ZONE);
}
- /**
- * @param p WGS84 lat/lon (ellipsoid GRS80) (in degree)
- * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
- */
- public EastNorth latlon2eastNorth(LatLon p) {
- double lt = Math.toRadians(p.lat());
- double lg = Math.toRadians(p.lon());
- if (lt >= cMinLatZonesRadian && lt <= cMaxLatZonesRadian && lg >= cMinLonZonesRadian && lg <= cMaxLonZonesRadian)
- return ConicProjection(lt, lg, layoutZone);
- return ConicProjection(lt, lg, 0);
+ public LambertCC9Zones(int layoutZone) {
+ updateParameters(layoutZone);
}
- /**
- *
- * @param lat latitude in grad
- * @param lon longitude in grad
- * @param nz Lambert CC zone number (from 1 to 9) - 1 !
- * @return EastNorth projected coordinates in meter
- */
- private EastNorth ConicProjection(double lat, double lon, int nz) {
- double R = c[nz]*Math.exp(-n[nz]*L(lat,e));
- double gamma = n[nz]*(lon-lambda0);
- double X = Xs + R*Math.sin(gamma);
- double Y = Ys[nz] + -R*Math.cos(gamma);
- return new EastNorth(X, Y);
+ public void updateParameters(int layoutZone) {
+ ellps = Ellipsoid.GRS80;
+ datum = GRS80Datum.INSTANCE;
+ this.layoutZone = layoutZone;
+ x_0 = 1700000;
+ y_0 = (layoutZone+1) * 1000000 + 200000;
+ lon_0 = 3;
+ double lat_0 = 42 + layoutZone;
+ double lat_1 = 41.25 + layoutZone;
+ double lat_2 = 42.75 + layoutZone;
+ if (proj == null) {
+ proj = new LambertConformalConic();
+ }
+ ((LambertConformalConic)proj).updateParameters2SP(ellps, lat_0, lat_1, lat_2);
}
- public LatLon eastNorth2latlon(EastNorth p) {
- return Geographic(p, layoutZone);
- }
-
- private LatLon Geographic(EastNorth ea, int nz) {
- double R = Math.sqrt(Math.pow(ea.getX()-Xs,2)+Math.pow(ea.getY()-Ys[nz], 2));
- double gamma = Math.atan((ea.getX()-Xs)/(Ys[nz]-ea.getY()));
- double lon = lambda0+gamma/n[nz];
- double latIso = (-1/n[nz])*Math.log(Math.abs(R/c[nz]));
- double lat = Ellipsoid.GRS80.latitude(latIso, e, epsilon);
- return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon));
- }
-
- @Override public String toString() {
+ @Override
+ public String toString() {
return tr("Lambert CC9 Zone (France)");
}
@@ -140,6 +77,7 @@
}
- public String toCode() {
- return "EPSG:"+(3942+layoutZone); //CC42 is EPSG:3942 (up to EPSG:3950 for CC50)
+ @Override
+ public Integer getEpsgCode() {
+ return 3942+layoutZone; //CC42 is EPSG:3942 (up to EPSG:3950 for CC50)
}
@@ -149,16 +87,10 @@
}
+ @Override
public String getCacheDirectoryName() {
return "lambert";
}
- /**
- * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
- */
- public double getDefaultZoomInPPD() {
- // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
- return 10.0;
- }
-
+ @Override
public Bounds getWorldBoundsLatLon()
{
@@ -203,15 +135,17 @@
}
+ @Override
public Collection getPreferences(JPanel p) {
Object prefcb = p.getComponent(2);
if(!(prefcb instanceof JComboBox))
return null;
- layoutZone = ((JComboBox)prefcb).getSelectedIndex();
+ int layoutZone = ((JComboBox)prefcb).getSelectedIndex();
return Collections.singleton(Integer.toString(layoutZone+1));
}
+ @Override
public void setPreferences(Collection args)
{
- layoutZone = DEFAULT_ZONE;
+ int layoutZone = DEFAULT_ZONE;
if (args != null) {
try {
@@ -226,4 +160,5 @@
} catch(NumberFormatException e) {}
}
+ updateParameters(layoutZone);
}
@@ -237,4 +172,5 @@
}
+ @Override
public Collection getPreferencesFromCode(String code)
{
@@ -246,5 +182,5 @@
if(zoneval >= 0 && zoneval <= 8)
return Collections.singleton(String.valueOf(zoneval+1));
- } catch(NumberFormatException e) {}
+ } catch(NumberFormatException ex) {}
}
return null;
Index: trunk/src/org/openstreetmap/josm/data/projection/LambertEST.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/LambertEST.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/LambertEST.java (revision 4285)
@@ -1,6 +1,3 @@
-//License: GPL. For details, see LICENSE file.
-//Thanks to Johan Montagnat and its geoconv java converter application
-//(http://www.i3s.unice.fr/~johan/gps/ , published under GPL license)
-//from which some code and constants have been reused here.
+// License: GPL. For details, see LICENSE file.
package org.openstreetmap.josm.data.projection;
@@ -8,86 +5,28 @@
import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
+import org.openstreetmap.josm.data.projection.proj.LambertConformalConic;
-public class LambertEST implements Projection {
+/**
+ * Estonian Coordinate System of 1997.
+ *
+ * Thanks to Johan Montagnat and its geoconv java converter application
+ * (http://www.i3s.unice.fr/~johan/gps/ , published under GPL license)
+ * from which some code and constants have been reused here.
+ */
+public class LambertEST extends AbstractProjection {
- public static final double ef = 500000; //Easting at false origin = 5000000 m
- public static final double nf = 6375000; //Northing at false origin = 6375000 m
- public static final double lat1 = Math.toRadians(59 + 1.0/3.0); //Latitude of 1st standard parallel = 59o20`0`` N
- public static final double lat2 = Math.toRadians( 58);//Latitude of 2nd standard parallel = 58o0`0`` N
- public static final double latf = Math.toRadians(57.517553930555555555555555555556);//'Latitude of false origin = 57o31`3,19415`` N
- public static final double lonf = Math.toRadians( 24.0);
- public static final double a = 6378137;
- public static final double ee = 0.081819191042815792;
- public static final double m1 = Math.cos(lat1) / (Math.sqrt(1 - ee *ee * Math.pow(Math.sin(lat1), 2)));
- public static final double m2 = Math.cos(lat2) / (Math.sqrt(1 - ee *ee * Math.pow(Math.sin(lat2), 2)));
- 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));
- 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));
- 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));
- public static final double n = (Math.log(m1) - Math.log(m2))
- / (Math.log(t1) - Math.log(t2));
- public static final double f = m1 / (n * Math.pow(t1, n));
- public static final double rf = a * f * Math.pow(tf, n);
-
- /**
- * precision in iterative schema
- */
- public static final double epsilon = 1e-11;
-
- /**
- * @param p WGS84 lat/lon (ellipsoid GRS80) (in degree)
- * @return eastnorth projection in Lambert Zone (ellipsoid GRS80)
- */
- public EastNorth latlon2eastNorth(LatLon p)
- {
-
- 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));
- double r = a * f * Math.pow(t, n);
- double theta = n * (Math.toRadians(p.lon()) - lonf);
-
- double x = ef + r * Math.sin(theta); //587446.7
- double y = nf + rf - r * Math.cos(theta); //6485401.6
-
- return new EastNorth(x,y);
- }
-
- public static double iterateAngle(double e, double t)
- {
- double a1 = 0.0;
- double a2 = 3.1415926535897931;
- double a = 1.5707963267948966;
- 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)));
- while (Math.abs(a-b) > epsilon)
- {
- a = a1 + ((a2 - a1) / 2.0);
- 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)));
- if (a1 == a2)
- return 0.0;
- if (b > a) {
- a1 = a;
- } else {
- a2 = a;
- }
- }
- return b;
- }
-
- public LatLon eastNorth2latlon(EastNorth p)
- {
- double r = Math.sqrt(Math.pow((p.getX() - ef), 2.0) + Math.pow((rf
- - p.getY() + nf), 2.0) ) * Math.signum(n);
- double T = Math.pow((r / (a * f)), (1.0/ n)) ;
- double theta = Math.atan((p.getX() - ef) / (rf - p.getY() + nf));
- double y = (theta / n + lonf) ;
- double x = iterateAngle(ee, T);
- return new LatLon(Math.toDegrees(x),Math.toDegrees(y));
+ public LambertEST() {
+ ellps = Ellipsoid.GRS80;
+ datum = GRS80Datum.INSTANCE;
+ lon_0 = 24;
+ double lat_0 = 57.517553930555555555555555555556;
+ double lat_1 = 59 + 1.0/3.0;
+ double lat_2 = 58;
+ x_0 = 500000;
+ y_0 = 6375000;
+ proj = new LambertConformalConic();
+ ((LambertConformalConic) proj).updateParameters2SP(ellps, lat_0, lat_1, lat_2);
}
@@ -97,6 +36,7 @@
}
- public String toCode() {
- return "EPSG:3301";
+ @Override
+ public Integer getEpsgCode() {
+ return 3301;
}
@@ -106,8 +46,10 @@
}
+ @Override
public String getCacheDirectoryName() {
return "lambertest";
}
+ @Override
public Bounds getWorldBoundsLatLon()
{
@@ -117,6 +59,3 @@
}
- public double getDefaultZoomInPPD() {
- return 10;
- }
}
Index: trunk/src/org/openstreetmap/josm/data/projection/Mercator.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/Mercator.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/Mercator.java (revision 4285)
@@ -5,10 +5,9 @@
import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.WGS84Datum;
/**
- * Implement Mercator Projection code, coded after documentation
- * from wikipedia.
+ * Mercator Projection
*
* The center of the mercator projection is always the 0 grad
@@ -20,26 +19,22 @@
* @author imi
*/
-public class Mercator implements Projection {
+public class Mercator extends AbstractProjection {
- final double radius = 6378137.0;
-
- public EastNorth latlon2eastNorth(LatLon p) {
- return new EastNorth(
- p.lon()*Math.PI/180*radius,
- Math.log(Math.tan(Math.PI/4+p.lat()*Math.PI/360))*radius);
+ public Mercator() {
+ ellps = Ellipsoid.WGS84;
+ datum = WGS84Datum.INSTANCE;
+ proj = new org.openstreetmap.josm.data.projection.proj.Mercator();
}
- public LatLon eastNorth2latlon(EastNorth p) {
- return new LatLon(
- Math.atan(Math.sinh(p.north()/radius))*180/Math.PI,
- p.east()/radius*180/Math.PI);
- }
-
- @Override public String toString() {
+ @Override
+ public String toString() {
return tr("Mercator");
}
- public String toCode() {
- return "EPSG:3857"; /* initially they used 3785 but that has been superseded, see http://www.epsg-registry.org/ */
+ @Override
+ public Integer getEpsgCode() {
+ /* initially they used 3785 but that has been superseded,
+ * see http://www.epsg-registry.org/ */
+ return 3857;
}
@@ -49,8 +44,10 @@
}
+ @Override
public String getCacheDirectoryName() {
return "mercator";
}
+ @Override
public Bounds getWorldBoundsLatLon()
{
@@ -60,7 +57,3 @@
}
- public double getDefaultZoomInPPD() {
- // This will set the scale bar to about 100 km
- return 1000.0;/*0.000158*/
- }
}
Index: trunk/src/org/openstreetmap/josm/data/projection/Puwg.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/Puwg.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/Puwg.java (revision 4285)
@@ -7,5 +7,4 @@
import java.awt.GridBagLayout;
import java.awt.event.ActionListener;
-import java.text.DecimalFormat;
import java.util.Collection;
import java.util.Collections;
@@ -16,6 +15,6 @@
import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
import org.openstreetmap.josm.tools.GBC;
@@ -26,9 +25,11 @@
* @author steelman
*/
-public class Puwg extends UTM implements Projection,ProjectionSubPrefs {
+public class Puwg extends AbstractProjection implements ProjectionSubPrefs {
+
public static final int DEFAULT_ZONE = 0;
- private int zone = DEFAULT_ZONE;
-
- static PuwgData[] Zones = new PuwgData[]{
+
+ private int zone;
+
+ static PuwgData[] Zones = new PuwgData[] {
new Epsg2180(),
new Epsg2176(),
@@ -38,34 +39,32 @@
};
- private static DecimalFormat decFormatter = new DecimalFormat("###0");
-
- @Override
- public EastNorth latlon2eastNorth(LatLon p) {
+ public Puwg() {
+ this(DEFAULT_ZONE);
+ }
+
+ public Puwg(int zone) {
+ ellps = Ellipsoid.GRS80;
+ proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
+ datum = GRS80Datum.INSTANCE;
+ updateParameters(zone);
+ }
+
+ public void updateParameters(int zone) {
+ this.zone = zone;
PuwgData z = Zones[zone];
- double easting = z.getPuwgFalseEasting();
- double northing = z.getPuwgFalseNorthing();
- double scale = z.getPuwgScaleFactor();
- double center = z.getPuwgCentralMeridian(); /* in radians */
- EastNorth a = mapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), center);
- return new EastNorth(a.east() * scale + easting, a.north() * scale + northing);
- }
-
- @Override
- public LatLon eastNorth2latlon(EastNorth p) {
- PuwgData z = Zones[zone];
- double easting = z.getPuwgFalseEasting();
- double northing = z.getPuwgFalseNorthing();
- double scale = z.getPuwgScaleFactor();
- double center = z.getPuwgCentralMeridian(); /* in radians */
- return mapXYToLatLon((p.east() - easting)/scale, (p.north() - northing)/scale, center);
- }
-
- @Override public String toString() {
+ x_0 = z.getPuwgFalseEasting();
+ y_0 = z.getPuwgFalseNorthing();
+ lon_0 = z.getPuwgCentralMeridianDeg();
+ k_0 = z.getPuwgScaleFactor();
+ }
+
+ @Override
+ public String toString() {
return tr("PUWG (Poland)");
}
@Override
- public String toCode() {
- return Zones[zone].toCode();
+ public Integer getEpsgCode() {
+ return Zones[zone].getEpsgCode();
}
@@ -83,18 +82,4 @@
public Bounds getWorldBoundsLatLon() {
return Zones[zone].getWorldBoundsLatLon();
- }
-
- @Override
- public double getDefaultZoomInPPD() {
- // This will set the scale bar to about 100 km
- return 0.009;
- }
-
- public String eastToString(EastNorth p) {
- return decFormatter.format(p.east());
- }
-
- public String northToString(EastNorth p) {
- return decFormatter.format(p.north());
}
@@ -135,9 +120,7 @@
@Override
- public Collection getPreferencesFromCode(String code)
- {
- for (PuwgData p : Puwg.Zones)
- {
- if(code.equals(p.toCode()))
+ public Collection getPreferencesFromCode(String code) {
+ for (PuwgData p : Puwg.Zones) {
+ if (code.equals(p.toCode()))
return Collections.singleton(code);
}
@@ -146,15 +129,13 @@
@Override
- public void setPreferences(Collection args)
- {
- zone = DEFAULT_ZONE;
- if(args != null)
- {
+ public void setPreferences(Collection args) {
+ int z = DEFAULT_ZONE;
+ if (args != null) {
try {
- for(String s : args)
- {
- for (int i=0; i < Puwg.Zones.length; ++i)
- if(s.equals(Zones[i].toCode())) {
- zone = i;
+ for (String s : args) {
+ for (int i=0; i < Zones.length; ++i)
+ if (s.equals(Zones[i].toCode())) {
+ z = i;
+ break;
}
break;
@@ -162,4 +143,5 @@
} catch (NullPointerException e) {}
}
+ updateParameters(z);
}
}
@@ -173,4 +155,5 @@
// Projection methods
+ Integer getEpsgCode();
String toCode();
String getCacheDirectoryName();
@@ -189,12 +172,20 @@
}
+ @Override
+ public Integer getEpsgCode() {
+ return 2180;
+ }
+
+ @Override
public String toCode() {
- return "EPSG:2180";
- }
-
+ return "EPSG:" + getEpsgCode();
+ }
+
+ @Override
public String getCacheDirectoryName() {
return "epsg2180";
}
+ @Override
public Bounds getWorldBoundsLatLon()
{
@@ -204,9 +195,9 @@
}
- public double getPuwgCentralMeridianDeg() { return Epsg2180CentralMeridian; }
- public double getPuwgCentralMeridian() { return Math.toRadians(Epsg2180CentralMeridian); }
- public double getPuwgFalseEasting() { return Epsg2180FalseEasting; }
- public double getPuwgFalseNorthing() { return Epsg2180FalseNorthing; }
- public double getPuwgScaleFactor() { return Epsg2180ScaleFactor; }
+ @Override public double getPuwgCentralMeridianDeg() { return Epsg2180CentralMeridian; }
+ @Override public double getPuwgCentralMeridian() { return Math.toRadians(Epsg2180CentralMeridian); }
+ @Override public double getPuwgFalseEasting() { return Epsg2180FalseEasting; }
+ @Override public double getPuwgFalseNorthing() { return Epsg2180FalseNorthing; }
+ @Override public double getPuwgScaleFactor() { return Epsg2180ScaleFactor; }
}
@@ -217,6 +208,6 @@
private static final double PuwgScaleFactor = 0.999923;
//final private double[] Puwg2000CentralMeridian = {15.0, 18.0, 21.0, 24.0};
- final private String[] Puwg2000Code = { "EPSG:2176", "EPSG:2177", "EPSG:2178", "EPSG:2179"};
- final private String[] Puwg2000CDName = { "epsg2176", "epsg2177", "epsg2178", "epsg2179"};
+ final private Integer[] Puwg2000Code = { 2176, 2177, 2178, 2179 };
+ final private String[] Puwg2000CDName = { "epsg2176", "epsg2177", "epsg2178", "epsg2179" };
@Override public String toString() {
@@ -224,12 +215,20 @@
}
+ @Override
+ public Integer getEpsgCode() {
+ return Puwg2000Code[getZoneIndex()];
+ }
+
+ @Override
public String toCode() {
- return Puwg2000Code[getZoneIndex()];
- }
-
+ return "EPSG:" + getEpsgCode();
+ }
+
+ @Override
public String getCacheDirectoryName() {
return Puwg2000CDName[getZoneIndex()];
}
+ @Override
public Bounds getWorldBoundsLatLon()
{
@@ -239,9 +238,9 @@
}
- public double getPuwgCentralMeridianDeg() { return getZone() * 3.0; }
- public double getPuwgCentralMeridian() { return Math.toRadians(getZone() * 3.0); }
- public double getPuwgFalseNorthing() { return PuwgFalseNorthing;}
- public double getPuwgFalseEasting() { return 1e6 * getZone() + PuwgFalseEasting; }
- public double getPuwgScaleFactor() { return PuwgScaleFactor; }
+ @Override public double getPuwgCentralMeridianDeg() { return getZone() * 3.0; }
+ @Override public double getPuwgCentralMeridian() { return Math.toRadians(getZone() * 3.0); }
+ @Override public double getPuwgFalseNorthing() { return PuwgFalseNorthing;}
+ @Override public double getPuwgFalseEasting() { return 1e6 * getZone() + PuwgFalseEasting; }
+ @Override public double getPuwgScaleFactor() { return PuwgScaleFactor; }
public abstract int getZone();
Index: trunk/src/org/openstreetmap/josm/data/projection/SwissGrid.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/SwissGrid.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/SwissGrid.java (revision 4285)
@@ -12,111 +12,32 @@
import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.ThreeParameterDatum;
+import org.openstreetmap.josm.data.projection.proj.SwissObliqueMercator;
import org.openstreetmap.josm.gui.widgets.HtmlPanel;
import org.openstreetmap.josm.tools.GBC;
/**
- * Projection for the SwissGrid CH1903 / L03, see http://de.wikipedia.org/wiki/Swiss_Grid.
+ * SwissGrid CH1903 / L03, see http://de.wikipedia.org/wiki/Swiss_Grid.
+ *
+ * Actually, what we have here, is CH1903+ (EPSG:2056), but without
+ * the additional false easting of 2000km and false northing 1000 km.
*
- * Calculations are based on formula from
- * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf
- *
- * August 2010 update to this formula (rigorous formulas)
- * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf
+ * To get to CH1903, a shift file is required. So currently, there are errors
+ * up to 1.6m (depending on the location).
+ *
+ * This projection does not have any parameters, it only implements
+ * ProjectionSubPrefs to show a warning that the grid file correction is not done.
*/
-public class SwissGrid implements Projection, ProjectionSubPrefs {
+public class SwissGrid extends AbstractProjection implements ProjectionSubPrefs {
- private static final double dX = 674.374;
- private static final double dY = 15.056;
- private static final double dZ = 405.346;
-
- private static final double phi0 = Math.toRadians(46.0 + 57.0 / 60 + 8.66 / 3600);
- private static final double lambda0 = Math.toRadians(7.0 + 26.0 / 60 + 22.50 / 3600);
- private static final double R = Ellipsoid.Bessel1841.a * Math.sqrt(1 - Ellipsoid.Bessel1841.e2) / (1 - (Ellipsoid.Bessel1841.e2 * Math.pow(Math.sin(phi0), 2)));
- private static final double alpha = Math.sqrt(1 + (Ellipsoid.Bessel1841.eb2 * Math.pow(Math.cos(phi0), 4)));
- private static final double b0 = Math.asin(Math.sin(phi0) / alpha);
- private static final double K = Math.log(Math.tan(Math.PI / 4 + b0 / 2)) - alpha
- * Math.log(Math.tan(Math.PI / 4 + phi0 / 2)) + alpha * Ellipsoid.Bessel1841.e / 2
- * Math.log((1 + Ellipsoid.Bessel1841.e * Math.sin(phi0)) / (1 - Ellipsoid.Bessel1841.e * Math.sin(phi0)));
-
- private static final double xTrans = 200000;
- private static final double yTrans = 600000;
-
- private static final double DELTA_PHI = 1e-11;
-
- private LatLon correctEllipoideGSR80toBressel1841(LatLon coord) {
- double[] XYZ = Ellipsoid.WGS84.latLon2Cart(coord);
- XYZ[0] -= dX;
- XYZ[1] -= dY;
- XYZ[2] -= dZ;
- return Ellipsoid.Bessel1841.cart2LatLon(XYZ);
- }
-
- private LatLon correctEllipoideBressel1841toGRS80(LatLon coord) {
- double[] XYZ = Ellipsoid.Bessel1841.latLon2Cart(coord);
- XYZ[0] += dX;
- XYZ[1] += dY;
- XYZ[2] += dZ;
- return Ellipsoid.WGS84.cart2LatLon(XYZ);
- }
-
- /**
- * @param wgs WGS84 lat/lon (ellipsoid GRS80) (in degree)
- * @return eastnorth projection in Swiss National Grid (ellipsoid Bessel)
- */
- @Override
- public EastNorth latlon2eastNorth(LatLon wgs) {
- LatLon coord = correctEllipoideGSR80toBressel1841(wgs);
- double phi = Math.toRadians(coord.lat());
- double lambda = Math.toRadians(coord.lon());
-
- double S = alpha * Math.log(Math.tan(Math.PI / 4 + phi / 2)) - alpha * Ellipsoid.Bessel1841.e / 2
- * Math.log((1 + Ellipsoid.Bessel1841.e * Math.sin(phi)) / (1 - Ellipsoid.Bessel1841.e * Math.sin(phi))) + K;
- double b = 2 * (Math.atan(Math.exp(S)) - Math.PI / 4);
- double l = alpha * (lambda - lambda0);
-
- double lb = Math.atan2(Math.sin(l), Math.sin(b0) * Math.tan(b) + Math.cos(b0) * Math.cos(l));
- double bb = Math.asin(Math.cos(b0) * Math.sin(b) - Math.sin(b0) * Math.cos(b) * Math.cos(l));
-
- double y = R * lb;
- double x = R / 2 * Math.log((1 + Math.sin(bb)) / (1 - Math.sin(bb)));
-
- return new EastNorth(y + yTrans, x + xTrans);
- }
-
- /**
- * @param xy SwissGrid east/north (in meters)
- * @return LatLon WGS84 (in degree)
- */
- @Override
- public LatLon eastNorth2latlon(EastNorth xy) {
- double x = xy.north() - xTrans;
- double y = xy.east() - yTrans;
-
- double lb = y / R;
- double bb = 2 * (Math.atan(Math.exp(x / R)) - Math.PI / 4);
-
- double b = Math.asin(Math.cos(b0) * Math.sin(bb) + Math.sin(b0) * Math.cos(bb) * Math.cos(lb));
- double l = Math.atan2(Math.sin(lb), Math.cos(b0) * Math.cos(lb) - Math.sin(b0) * Math.tan(bb));
-
- double lambda = lambda0 + l / alpha;
- double phi = b;
- double S = 0;
-
- double prevPhi = -1000;
- int iteration = 0;
- // iteration to finds S and phi
- while (Math.abs(phi - prevPhi) > DELTA_PHI) {
- if (++iteration > 30)
- throw new RuntimeException("Two many iterations");
- prevPhi = phi;
- S = 1 / alpha * (Math.log(Math.tan(Math.PI / 4 + b / 2)) - K) + Ellipsoid.Bessel1841.e
- * Math.log(Math.tan(Math.PI / 4 + Math.asin(Ellipsoid.Bessel1841.e * Math.sin(phi)) / 2));
- phi = 2 * Math.atan(Math.exp(S)) - Math.PI / 2;
- }
-
- LatLon coord = correctEllipoideBressel1841toGRS80(new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda)));
- return coord;
+ public SwissGrid() {
+ ellps = Ellipsoid.Bessel1841;
+ datum = new ThreeParameterDatum("SwissGrid Datum", null, ellps, 674.374, 15.056, 405.346);
+ x_0 = 600000;
+ y_0 = 200000;
+ lon_0 = 7.0 + 26.0/60 + 22.50/3600;
+ double lat_0 = 46.0 + 57.0/60 + 8.66/3600;
+ proj = new SwissObliqueMercator(ellps, lat_0);
}
@@ -127,6 +48,6 @@
@Override
- public String toCode() {
- return "EPSG:21781";
+ public Integer getEpsgCode() {
+ return 21781;
}
@@ -144,10 +65,4 @@
public Bounds getWorldBoundsLatLon() {
return new Bounds(new LatLon(45.7, 5.7), new LatLon(47.9, 10.6));
- }
-
- @Override
- public double getDefaultZoomInPPD() {
- // This will set the scale bar to about 100 m
- return 1.01;
}
Index: trunk/src/org/openstreetmap/josm/data/projection/TransverseMercator.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/TransverseMercator.java (revision 4284)
+++ (revision )
@@ -1,328 +1,0 @@
-// License: GPL. For details, see LICENSE file.
-package org.openstreetmap.josm.data.projection;
-
-import org.openstreetmap.josm.data.coor.EastNorth;
-import org.openstreetmap.josm.data.coor.LatLon;
-
-/**
- * This is a base class to do projections based on Transverse Mercator projection.
- *
- * @author Dirk Stöcker
- * code based on JavaScript from Chuck Taylor
- *
- * NOTE: Uses polygon approximation to translate to WGS84.
- */
-public abstract class TransverseMercator implements Projection {
-
- private final static double UTMScaleFactor = 0.9996;
-
- private double UTMCentralMeridianRad = 0;
- private double offsetEastMeters = 500000;
- private double offsetNorthMeters = 0;
-
-
- protected void setProjectionParameters(double centralMeridianDegress, double offsetEast, double offsetNorth)
- {
- UTMCentralMeridianRad = Math.toRadians(centralMeridianDegress);
- offsetEastMeters = offsetEast;
- offsetNorthMeters = offsetNorth;
- }
-
- /*
- * ArcLengthOfMeridian
- *
- * Computes the ellipsoidal distance from the equator to a point at a
- * given latitude.
- *
- * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
- * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
- *
- * Inputs:
- * phi - Latitude of the point, in radians.
- *
- * Globals:
- * Ellipsoid.GRS80.a - Ellipsoid model major axis.
- * Ellipsoid.GRS80.b - Ellipsoid model minor axis.
- *
- * Returns:
- * The ellipsoidal distance of the point from the equator, in meters.
- *
- */
- private double ArcLengthOfMeridian(double phi)
- {
- /* Precalculate n */
- double n = (Ellipsoid.GRS80.a - Ellipsoid.GRS80.b) / (Ellipsoid.GRS80.a + Ellipsoid.GRS80.b);
-
- /* Precalculate alpha */
- double alpha = ((Ellipsoid.GRS80.a + Ellipsoid.GRS80.b) / 2.0)
- * (1.0 + (Math.pow (n, 2.0) / 4.0) + (Math.pow (n, 4.0) / 64.0));
-
- /* Precalculate beta */
- double beta = (-3.0 * n / 2.0) + (9.0 * Math.pow (n, 3.0) / 16.0)
- + (-3.0 * Math.pow (n, 5.0) / 32.0);
-
- /* Precalculate gamma */
- double gamma = (15.0 * Math.pow (n, 2.0) / 16.0)
- + (-15.0 * Math.pow (n, 4.0) / 32.0);
-
- /* Precalculate delta */
- double delta = (-35.0 * Math.pow (n, 3.0) / 48.0)
- + (105.0 * Math.pow (n, 5.0) / 256.0);
-
- /* Precalculate epsilon */
- double epsilon = (315.0 * Math.pow (n, 4.0) / 512.0);
-
- /* Now calculate the sum of the series and return */
- return alpha
- * (phi + (beta * Math.sin (2.0 * phi))
- + (gamma * Math.sin (4.0 * phi))
- + (delta * Math.sin (6.0 * phi))
- + (epsilon * Math.sin (8.0 * phi)));
- }
-
- /*
- * FootpointLatitude
- *
- * Computes the footpoint latitude for use in converting transverse
- * Mercator coordinates to ellipsoidal coordinates.
- *
- * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
- * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
- *
- * Inputs:
- * y - The UTM northing coordinate, in meters.
- *
- * Returns:
- * The footpoint latitude, in radians.
- *
- */
- private double FootpointLatitude(double y)
- {
- /* Precalculate n (Eq. 10.18) */
- double n = (Ellipsoid.GRS80.a - Ellipsoid.GRS80.b) / (Ellipsoid.GRS80.a + Ellipsoid.GRS80.b);
-
- /* Precalculate alpha_ (Eq. 10.22) */
- /* (Same as alpha in Eq. 10.17) */
- double alpha_ = ((Ellipsoid.GRS80.a + Ellipsoid.GRS80.b) / 2.0)
- * (1 + (Math.pow (n, 2.0) / 4) + (Math.pow (n, 4.0) / 64));
-
- /* Precalculate y_ (Eq. 10.23) */
- double y_ = y / alpha_;
-
- /* Precalculate beta_ (Eq. 10.22) */
- double beta_ = (3.0 * n / 2.0) + (-27.0 * Math.pow (n, 3.0) / 32.0)
- + (269.0 * Math.pow (n, 5.0) / 512.0);
-
- /* Precalculate gamma_ (Eq. 10.22) */
- double gamma_ = (21.0 * Math.pow (n, 2.0) / 16.0)
- + (-55.0 * Math.pow (n, 4.0) / 32.0);
-
- /* Precalculate delta_ (Eq. 10.22) */
- double delta_ = (151.0 * Math.pow (n, 3.0) / 96.0)
- + (-417.0 * Math.pow (n, 5.0) / 128.0);
-
- /* Precalculate epsilon_ (Eq. 10.22) */
- double epsilon_ = (1097.0 * Math.pow (n, 4.0) / 512.0);
-
- /* Now calculate the sum of the series (Eq. 10.21) */
- return y_ + (beta_ * Math.sin (2.0 * y_))
- + (gamma_ * Math.sin (4.0 * y_))
- + (delta_ * Math.sin (6.0 * y_))
- + (epsilon_ * Math.sin (8.0 * y_));
- }
-
- /*
- * MapLatLonToXY
- *
- * Converts a latitude/longitude pair to x and y coordinates in the
- * Transverse Mercator projection. Note that Transverse Mercator is not
- * the same as UTM; a scale factor is required to convert between them.
- *
- * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
- * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
- *
- * Inputs:
- * phi - Latitude of the point, in radians.
- * lambda - Longitude of the point, in radians.
- * lambda0 - Longitude of the central meridian to be used, in radians.
- *
- * Outputs:
- * xy - A 2-element array containing the x and y coordinates
- * of the computed point.
- *
- * Returns:
- * The function does not return a value.
- *
- */
- public EastNorth mapLatLonToXY(double phi, double lambda, double lambda0)
- {
- /* Precalculate ep2 */
- double ep2 = (Math.pow (Ellipsoid.GRS80.a, 2.0) - Math.pow (Ellipsoid.GRS80.b, 2.0)) / Math.pow (Ellipsoid.GRS80.b, 2.0);
-
- /* Precalculate nu2 */
- double nu2 = ep2 * Math.pow (Math.cos (phi), 2.0);
-
- /* Precalculate N */
- double N = Math.pow (Ellipsoid.GRS80.a, 2.0) / (Ellipsoid.GRS80.b * Math.sqrt (1 + nu2));
-
- /* Precalculate t */
- double t = Math.tan (phi);
- double t2 = t * t;
-
- /* Precalculate l */
- double l = lambda - lambda0;
-
- /* Precalculate coefficients for l**n in the equations below
- so a normal human being can read the expressions for easting
- and northing
- -- l**1 and l**2 have coefficients of 1.0 */
- double l3coef = 1.0 - t2 + nu2;
-
- double l4coef = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2);
-
- double l5coef = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2
- - 58.0 * t2 * nu2;
-
- double l6coef = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2
- - 330.0 * t2 * nu2;
-
- double l7coef = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2);
-
- double l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
-
- return new EastNorth(
- /* Calculate easting (x) */
- N * Math.cos (phi) * l
- + (N / 6.0 * Math.pow (Math.cos (phi), 3.0) * l3coef * Math.pow (l, 3.0))
- + (N / 120.0 * Math.pow (Math.cos (phi), 5.0) * l5coef * Math.pow (l, 5.0))
- + (N / 5040.0 * Math.pow (Math.cos (phi), 7.0) * l7coef * Math.pow (l, 7.0)),
- /* Calculate northing (y) */
- ArcLengthOfMeridian (phi)
- + (t / 2.0 * N * Math.pow (Math.cos (phi), 2.0) * Math.pow (l, 2.0))
- + (t / 24.0 * N * Math.pow (Math.cos (phi), 4.0) * l4coef * Math.pow (l, 4.0))
- + (t / 720.0 * N * Math.pow (Math.cos (phi), 6.0) * l6coef * Math.pow (l, 6.0))
- + (t / 40320.0 * N * Math.pow (Math.cos (phi), 8.0) * l8coef * Math.pow (l, 8.0)));
- }
-
- /*
- * MapXYToLatLon
- *
- * Converts x and y coordinates in the Transverse Mercator projection to
- * a latitude/longitude pair. Note that Transverse Mercator is not
- * the same as UTM; a scale factor is required to convert between them.
- *
- * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
- * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
- *
- * Inputs:
- * x - The easting of the point, in meters.
- * y - The northing of the point, in meters.
- * lambda0 - Longitude of the central meridian to be used, in radians.
- *
- * Outputs:
- * philambda - A 2-element containing the latitude and longitude
- * in radians.
- *
- * Returns:
- * The function does not return a value.
- *
- * Remarks:
- * The local variables Nf, nuf2, tf, and tf2 serve the same purpose as
- * N, nu2, t, and t2 in MapLatLonToXY, but they are computed with respect
- * to the footpoint latitude phif.
- *
- * x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
- * to optimize computations.
- *
- */
- public LatLon mapXYToLatLon(double x, double y, double lambda0)
- {
- /* Get the value of phif, the footpoint latitude. */
- double phif = FootpointLatitude (y);
-
- /* Precalculate ep2 */
- double ep2 = (Math.pow (Ellipsoid.GRS80.a, 2.0) - Math.pow (Ellipsoid.GRS80.b, 2.0))
- / Math.pow (Ellipsoid.GRS80.b, 2.0);
-
- /* Precalculate cos (phif) */
- double cf = Math.cos (phif);
-
- /* Precalculate nuf2 */
- double nuf2 = ep2 * Math.pow (cf, 2.0);
-
- /* Precalculate Nf and initialize Nfpow */
- double Nf = Math.pow (Ellipsoid.GRS80.a, 2.0) / (Ellipsoid.GRS80.b * Math.sqrt (1 + nuf2));
- double Nfpow = Nf;
-
- /* Precalculate tf */
- double tf = Math.tan (phif);
- double tf2 = tf * tf;
- double tf4 = tf2 * tf2;
-
- /* Precalculate fractional coefficients for x**n in the equations
- below to simplify the expressions for latitude and longitude. */
- double x1frac = 1.0 / (Nfpow * cf);
-
- Nfpow *= Nf; /* now equals Nf**2) */
- double x2frac = tf / (2.0 * Nfpow);
-
- Nfpow *= Nf; /* now equals Nf**3) */
- double x3frac = 1.0 / (6.0 * Nfpow * cf);
-
- Nfpow *= Nf; /* now equals Nf**4) */
- double x4frac = tf / (24.0 * Nfpow);
-
- Nfpow *= Nf; /* now equals Nf**5) */
- double x5frac = 1.0 / (120.0 * Nfpow * cf);
-
- Nfpow *= Nf; /* now equals Nf**6) */
- double x6frac = tf / (720.0 * Nfpow);
-
- Nfpow *= Nf; /* now equals Nf**7) */
- double x7frac = 1.0 / (5040.0 * Nfpow * cf);
-
- Nfpow *= Nf; /* now equals Nf**8) */
- double x8frac = tf / (40320.0 * Nfpow);
-
- /* Precalculate polynomial coefficients for x**n.
- -- x**1 does not have a polynomial coefficient. */
- double x2poly = -1.0 - nuf2;
- double x3poly = -1.0 - 2 * tf2 - nuf2;
- double x4poly = 5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 - 3.0 * (nuf2 *nuf2) - 9.0 * tf2 * (nuf2 * nuf2);
- double x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2;
- double x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 + 162.0 * tf2 * nuf2;
- double x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2);
- double x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2);
-
- return new LatLon(
- /* Calculate latitude */
- Math.toDegrees(
- phif + x2frac * x2poly * (x * x)
- + x4frac * x4poly * Math.pow (x, 4.0)
- + x6frac * x6poly * Math.pow (x, 6.0)
- + x8frac * x8poly * Math.pow (x, 8.0)),
- Math.toDegrees(
- /* Calculate longitude */
- lambda0 + x1frac * x
- + x3frac * x3poly * Math.pow (x, 3.0)
- + x5frac * x5poly * Math.pow (x, 5.0)
- + x7frac * x7poly * Math.pow (x, 7.0)));
- }
-
- @Override
- public EastNorth latlon2eastNorth(LatLon p) {
- EastNorth a = mapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), UTMCentralMeridianRad);
- return new EastNorth(a.east() * UTMScaleFactor + offsetEastMeters, a.north() * UTMScaleFactor + offsetNorthMeters);
- }
-
- @Override
- public LatLon eastNorth2latlon(EastNorth p) {
- return mapXYToLatLon((p.east() - offsetEastMeters)/UTMScaleFactor, (p.north() - offsetNorthMeters)/UTMScaleFactor, UTMCentralMeridianRad);
- }
-
- @Override
- public double getDefaultZoomInPPD() {
- // this will set the map scaler to about 1000 m
- return 10;
- }
-}
Index: trunk/src/org/openstreetmap/josm/data/projection/TransverseMercatorLV.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/TransverseMercatorLV.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/TransverseMercatorLV.java (revision 4285)
@@ -6,4 +6,5 @@
import org.openstreetmap.josm.data.Bounds;
import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
/**
@@ -13,22 +14,24 @@
* @author Viesturs Zarins
*/
-public class TransverseMercatorLV extends TransverseMercator {
+public class TransverseMercatorLV extends AbstractProjection {
- public TransverseMercatorLV()
- {
- setProjectionParameters(24, 500000, -6000000);
+ public TransverseMercatorLV() {
+ ellps = Ellipsoid.GRS80;
+ proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
+ datum = GRS80Datum.INSTANCE;
+ lon_0 = 24;
+ x_0 = 500000;
+ y_0 = -6000000;
+ k_0 = 0.9996;
}
-
- @Override public String toString() {
+
+ @Override
+ public String toString() {
return tr("LKS-92 (Latvia TM)");
}
- private int epsgCode() {
+ @Override
+ public Integer getEpsgCode() {
return 3059;
- }
-
- @Override
- public String toCode() {
- return "EPSG:"+ epsgCode();
}
@@ -38,6 +41,7 @@
}
+ @Override
public String getCacheDirectoryName() {
- return "epsg"+ epsgCode();
+ return "epsg"+ getEpsgCode();
}
Index: trunk/src/org/openstreetmap/josm/data/projection/UTM.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/UTM.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/UTM.java (revision 4285)
@@ -19,4 +19,5 @@
import org.openstreetmap.josm.data.Bounds;
import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
import org.openstreetmap.josm.tools.GBC;
@@ -25,19 +26,39 @@
* @author Dirk Stöcker
* code based on JavaScript from Chuck Taylor
+ *
*/
-public class UTM extends TransverseMercator implements ProjectionSubPrefs {
+public class UTM extends AbstractProjection implements ProjectionSubPrefs {
private static final int DEFAULT_ZONE = 30;
- private int zone = DEFAULT_ZONE;
-
- public enum Hemisphere {North, South}
+ private int zone;
+
+ public enum Hemisphere { North, South }
private static final Hemisphere DEFAULT_HEMISPHERE = Hemisphere.North;
- private Hemisphere hemisphere = DEFAULT_HEMISPHERE;
-
- private boolean offset = false;
-
- public UTM()
- {
- updateParameters();
+ private Hemisphere hemisphere;
+
+ /**
+ * Applies an additional false easting of 3000000 m if true.
+ */
+ private boolean offset;
+
+ public UTM() {
+ this(DEFAULT_ZONE, DEFAULT_HEMISPHERE, false);
+ }
+
+ public UTM(int zone, Hemisphere hemisphere, boolean offset) {
+ ellps = Ellipsoid.GRS80;
+ proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
+ datum = GRS80Datum.INSTANCE;
+ updateParameters(zone, hemisphere, offset);
+ }
+
+ public void updateParameters(int zone, Hemisphere hemisphere, boolean offset) {
+ this.zone = zone;
+ this.hemisphere = hemisphere;
+ this.offset = offset;
+ x_0 = 500000 + (offset ? 3000000 : 0);
+ y_0 = hemisphere == Hemisphere.North ? 0 : 10000000;
+ lon_0 = getUtmCentralMeridianDeg(zone);
+ k_0 = 0.9996;
}
@@ -56,39 +77,21 @@
*
*/
- private double UTMCentralMeridianDeg(int zone)
+ private double getUtmCentralMeridianDeg(int zone)
{
return -183.0 + (zone * 6.0);
}
- @Override public String toString() {
+ public int getzone() {
+ return zone;
+ }
+
+ @Override
+ public String toString() {
return tr("UTM");
}
- private void updateParameters() {
- setProjectionParameters(this.UTMCentralMeridianDeg(getzone()), getEastOffset(), getNorthOffset());
- }
-
- public int getzone()
- {
- return zone;
- }
-
- private double getEastOffset() {
- return 500000 + (offset?3000000:0);
- }
-
- private double getNorthOffset() {
- if (hemisphere == Hemisphere.North)
- return 0;
- else
- return 10000000;
- }
-
- private int epsgCode() {
+ @Override
+ public Integer getEpsgCode() {
return ((offset?325800:32600) + getzone() + (hemisphere == Hemisphere.South?100:0));
- }
-
- public String toCode() {
- return "EPSG:"+ epsgCode();
}
@@ -98,18 +101,20 @@
}
+ @Override
public String getCacheDirectoryName() {
- return "epsg"+ epsgCode();
- }
-
+ return "epsg"+ getEpsgCode();
+ }
+
+ @Override
public Bounds getWorldBoundsLatLon()
{
if (hemisphere == Hemisphere.North)
return new Bounds(
- new LatLon(-5.0, UTMCentralMeridianDeg(getzone())-5.0),
- new LatLon(85.0, UTMCentralMeridianDeg(getzone())+5.0));
+ new LatLon(-5.0, getUtmCentralMeridianDeg(getzone())-5.0),
+ new LatLon(85.0, getUtmCentralMeridianDeg(getzone())+5.0));
else
return new Bounds(
- new LatLon(-85.0, UTMCentralMeridianDeg(getzone())-5.0),
- new LatLon(5.0, UTMCentralMeridianDeg(getzone())+5.0));
+ new LatLon(-85.0, getUtmCentralMeridianDeg(getzone())-5.0),
+ new LatLon(5.0, getUtmCentralMeridianDeg(getzone())+5.0));
}
@@ -173,4 +178,5 @@
}
+ @Override
public Collection getPreferences(JPanel p) {
int zone = DEFAULT_ZONE;
@@ -199,9 +205,9 @@
}
- public void setPreferences(Collection args)
- {
- zone = DEFAULT_ZONE;
- hemisphere = DEFAULT_HEMISPHERE;
- offset = false;
+ @Override
+ public void setPreferences(Collection args) {
+ int zone = DEFAULT_ZONE;
+ Hemisphere hemisphere = DEFAULT_HEMISPHERE;
+ boolean offset = false;
if(args != null)
@@ -223,7 +229,8 @@
}
}
- updateParameters();
- }
-
+ updateParameters(zone, hemisphere, offset);
+ }
+
+ @Override
public String[] allCodes() {
ArrayList projections = new ArrayList(60*4);
@@ -236,9 +243,9 @@
}
return projections.toArray(new String[0]);
-
- }
-
- public Collection getPreferencesFromCode(String code)
- {
+ }
+
+ @Override
+ public Collection getPreferencesFromCode(String code) {
+
boolean offset = code.startsWith("EPSG:3258") || code.startsWith("EPSG:3259");
Index: trunk/src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java (revision 4284)
+++ trunk/src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java (revision 4285)
@@ -2,8 +2,4 @@
package org.openstreetmap.josm.data.projection;
-/**
- * This class implements all projections for French departements in the Caribbean Sea and
- * Indian Ocean using the UTM transvers Mercator projection and specific geodesic settings (7 parameters transformation algorithm).
- */
import static org.openstreetmap.josm.tools.I18n.tr;
@@ -18,44 +14,48 @@
import org.openstreetmap.josm.data.Bounds;
-import org.openstreetmap.josm.data.coor.EastNorth;
import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.datum.Datum;
+import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
+import org.openstreetmap.josm.data.projection.datum.SevenParameterDatum;
+import org.openstreetmap.josm.data.projection.datum.ThreeParameterDatum;
import org.openstreetmap.josm.tools.GBC;
-public class UTM_France_DOM implements Projection, ProjectionSubPrefs {
+/**
+ * This class implements all projections for French departements in the Caribbean Sea and
+ * Indian Ocean using the UTM transvers Mercator projection and specific geodesic settings (7 parameters transformation algorithm).
+ *
+ */
+public class UTM_France_DOM extends AbstractProjection implements ProjectionSubPrefs {
- private static String FortMarigotName = tr("Guadeloupe Fort-Marigot 1949");
- private static String SainteAnneName = tr("Guadeloupe Ste-Anne 1948");
- private static String MartiniqueName = tr("Martinique Fort Desaix 1952");
- private static String Reunion92Name = tr("Reunion RGR92");
- private static String Guyane92Name = tr("Guyane RGFG95");
- public static String[] utmGeodesicsNames = { FortMarigotName, SainteAnneName, MartiniqueName, Reunion92Name, Guyane92Name};
+ private final static String FortMarigotName = tr("Guadeloupe Fort-Marigot 1949");
+ private final static String SainteAnneName = tr("Guadeloupe Ste-Anne 1948");
+ private final static String MartiniqueName = tr("Martinique Fort Desaix 1952");
+ private final static String Reunion92Name = tr("Reunion RGR92");
+ private final static String Guyane92Name = tr("Guyane RGFG95");
+ private final static String[] utmGeodesicsNames = { FortMarigotName, SainteAnneName, MartiniqueName, Reunion92Name, Guyane92Name};
- private Bounds FortMarigotBounds = new Bounds( new LatLon(17.6,-63.25), new LatLon(18.5,-62.5));
- private Bounds SainteAnneBounds = new Bounds( new LatLon(15.8,-61.9), new LatLon(16.6,-60.9));
- private Bounds MartiniqueBounds = new Bounds( new LatLon(14.25,-61.25), new LatLon(15.025,-60.725));
- private Bounds ReunionBounds = new Bounds( new LatLon(-25.92,37.58), new LatLon(-10.6, 58.27));
- private Bounds GuyaneBounds = new Bounds( new LatLon(2.16 , -54.0), new LatLon(9.06 , -49.62));
- private Bounds[] utmBounds = { FortMarigotBounds, SainteAnneBounds, MartiniqueBounds, ReunionBounds, GuyaneBounds};
+ private final static Bounds FortMarigotBounds = new Bounds( new LatLon(17.6,-63.25), new LatLon(18.5,-62.5));
+ private final static Bounds SainteAnneBounds = new Bounds( new LatLon(15.8,-61.9), new LatLon(16.6,-60.9));
+ private final static Bounds MartiniqueBounds = new Bounds( new LatLon(14.25,-61.25), new LatLon(15.025,-60.725));
+ private final static Bounds ReunionBounds = new Bounds( new LatLon(-25.92,37.58), new LatLon(-10.6, 58.27));
+ private final static Bounds GuyaneBounds = new Bounds( new LatLon(2.16 , -54.0), new LatLon(9.06 , -49.62));
+ private final static Bounds[] utmBounds = { FortMarigotBounds, SainteAnneBounds, MartiniqueBounds, ReunionBounds, GuyaneBounds };
- private String FortMarigotEPSG = "EPSG::2969";
- private String SainteAnneEPSG = "EPSG::2970";
- private String MartiniqueEPSG = "EPSG::2973";
- private String ReunionEPSG = "EPSG::2975";
- private String GuyaneEPSG = "EPSG::2972";
- private String[] utmEPSGs = { FortMarigotEPSG, SainteAnneEPSG, MartiniqueEPSG, ReunionEPSG, GuyaneEPSG};
+ private final static Integer FortMarigotEPSG = 2969;
+ private final static Integer SainteAnneEPSG = 2970;
+ private final static Integer MartiniqueEPSG = 2973;
+ private final static Integer ReunionEPSG = 2975;
+ private final static Integer GuyaneEPSG = 2972;
+ private final static Integer[] utmEPSGs = { FortMarigotEPSG, SainteAnneEPSG, MartiniqueEPSG, ReunionEPSG, GuyaneEPSG };
- /**
- * false east in meters (constant)
- */
- private static final double Xs = 500000.0;
- /**
- * false north in meters (0 in northern hemisphere, 10000000 in southern
- * hemisphere)
- */
- private static double Ys = 0;
- /**
- * origin meridian longitude
- */
- protected double lg0;
+ private final static Datum FortMarigotDatum = new ThreeParameterDatum("FortMarigot Datum", null, Ellipsoid.hayford, 136.596, 248.148, -429.789);
+ private final static Datum SainteAnneDatum = new SevenParameterDatum("SainteAnne Datum", null, Ellipsoid.hayford, -472.29, -5.63, -304.12, 0.4362, -0.8374, 0.2563, 1.8984);
+ private final static Datum MartiniqueDatum = new SevenParameterDatum("Martinique Datum", null, Ellipsoid.hayford, 126.926, 547.939, 130.409, -2.78670, 5.16124, -0.85844, 13.82265);
+ private final static Datum ReunionDatum = GRS80Datum.INSTANCE;
+ private final static Datum GuyaneDatum = GRS80Datum.INSTANCE;
+ private final static Datum[] utmDatums = { FortMarigotDatum, SainteAnneDatum, MartiniqueDatum, ReunionDatum, GuyaneDatum };
+
+ private final static int[] utmZones = { 20, 20, 20, 40, 22 };
+
/**
* UTM zone (from 1 to 60)
@@ -69,330 +69,45 @@
public static final int DEFAULT_GEODESIC = 0;
- public static int currentGeodesic = DEFAULT_GEODESIC;
+ public int currentGeodesic;
- /**
- * 7 parameters transformation
- */
- private static double tx = 0.0;
- private static double ty = 0.0;
- private static double tz = 0.0;
- private static double rx = 0;
- private static double ry = 0;
- private static double rz = 0;
- private static double scaleDiff = 0;
- /**
- * precision in iterative schema
- */
- public static final double epsilon = 1e-11;
- private void refresh7ParametersTranslation() {
- if (currentGeodesic == 0) { // UTM_20N_Guadeloupe_Fort_Marigot
- set7ParametersTranslation(new double[]{136.596, 248.148, -429.789},
- new double[]{0, 0, 0},
- 0,
- true, 20);
- } else if (currentGeodesic == 1) { // UTM_20N_Guadeloupe_Ste_Anne
- set7ParametersTranslation(new double[]{-472.29, -5.63, -304.12},
- new double[]{0.4362, -0.8374, 0.2563},
- 1.8984E-6,
- true, 20);
- } else if (currentGeodesic == 2) { // UTM_20N_Martinique_Fort_Desaix
- set7ParametersTranslation(new double[]{126.926, 547.939, 130.409},
- new double[]{-2.78670, 5.16124, -0.85844},
- 13.82265E-6
- , true, 20);
- } else if (currentGeodesic == 3) { // UTM_40S_Reunion_RGR92 (translation only required for re-projections from Gauss-Laborde)
- set7ParametersTranslation(new double[]{789.524, -626.486, -89.904},
- new double[]{0.6006, 76.7946, -10.5788},
- -32.3241E-6
- , false, 40);
- } else if (currentGeodesic == 4) { // UTM_22N_Guyane_RGFG95 (translation only required for re-projections from CSG67)
- set7ParametersTranslation(new double[]{-193.066 , 236.993, 105.447},
- new double[]{0.4814, -0.8074, 0.1276},
- 1.5649E-6
- , true, 22);
- }
+ public UTM_France_DOM() {
+ updateParameters(DEFAULT_GEODESIC);
+ }
+
+ public void updateParameters(int currentGeodesic) {
+ this.currentGeodesic = currentGeodesic;
+ datum = utmDatums[currentGeodesic];
+ ellps = datum.getEllipsoid();
+ proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
+ isNorth = currentGeodesic != 3;
+ zone = utmZones[currentGeodesic];
+ x_0 = 500000;
+ y_0 = isNorth ? 0.0 : 10000000.0;
+ lon_0 = 6 * zone - 183;
+ k_0 = 0.9996;
+ }
+
+ public int getCurrentGeodesic() {
+ return currentGeodesic;
}
- private void set7ParametersTranslation(double[] translation, double[] rotation, double scalediff, boolean north, int utmZone) {
- tx = translation[0];
- ty = translation[1];
- tz = translation[2];
- rx = rotation[0]/206264.806247096355; // seconds to radian
- ry = rotation[1]/206264.806247096355;
- rz = rotation[2]/206264.806247096355;
- scaleDiff = scalediff;
- isNorth = north;
- Ys = isNorth ? 0.0 : 10000000.0;
- zone = utmZone;
+ @Override
+ public String toString() {
+ return tr("UTM France (DOM)");
}
- public EastNorth latlon2eastNorth(LatLon p) {
- if (currentGeodesic < 3 ) {
- // translate ellipsoid GRS80 (WGS83) => reference ellipsoid geographic
- LatLon geo = GRS802Hayford(p);
- // reference ellipsoid geographic => UTM projection
- return MTProjection(geo, Ellipsoid.hayford.a, Ellipsoid.hayford.e);
- } else { // UTM_40S_Reunion_RGR92 or UTM_22N_Guyane_RGFG95
- LatLon geo = new LatLon(Math.toRadians(p.lat()), Math.toRadians(p.lon()));
- return MTProjection(geo, Ellipsoid.GRS80.a, Ellipsoid.GRS80.e);
- }
- }
-
- /**
- * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to UTM
- * geographic, (ellipsoid Hayford)
- */
- private LatLon GRS802Hayford(LatLon wgs) {
- double lat = Math.toRadians(wgs.lat()); // degree to radian
- double lon = Math.toRadians(wgs.lon());
- // WGS84 geographic => WGS84 cartesian
- double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
- double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
- double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
- double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
- // translation
- double coord[] = invSevenParametersTransformation(X, Y, Z);
- // UTM cartesian => UTM geographic
- return Geographic(coord[0], coord[1], coord[2], Ellipsoid.hayford);
- }
-
- /**
- * initializes from cartesian coordinates
- *
- * @param X
- * 1st coordinate in meters
- * @param Y
- * 2nd coordinate in meters
- * @param Z
- * 3rd coordinate in meters
- * @param ell
- * reference ellipsoid
- */
- private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
- double norm = Math.sqrt(X * X + Y * Y);
- double lg = 2.0 * Math.atan(Y / (X + norm));
- double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
- double delta = 1.0;
- while (delta > epsilon) {
- double s2 = Math.sin(lt);
- s2 *= s2;
- double l = Math.atan((Z / norm)
- / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
- delta = Math.abs(l - lt);
- lt = l;
- }
- double s2 = Math.sin(lt);
- s2 *= s2;
- // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
- return new LatLon(lt, lg);
- }
-
- /**
- * initalizes from geographic coordinates
- *
- * @param coord geographic coordinates triplet
- * @param a reference ellipsoid long axis
- * @param e reference ellipsoid excentricity
- */
- private EastNorth MTProjection(LatLon coord, double a, double e) {
- double n = 0.9996 * a;
- Ys = (coord.lat() >= 0.0) ? 0.0 : 10000000.0;
- double r6d = Math.PI / 30.0;
- //zone = (int) Math.floor((coord.lon() + Math.PI) / r6d) + 1;
- lg0 = r6d * (zone - 0.5) - Math.PI;
- double e2 = e * e;
- double e4 = e2 * e2;
- double e6 = e4 * e2;
- double e8 = e4 * e4;
- double C[] = {
- 1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
- e2/8.0 - e4/96.0 - 9.0*e6/1024.0 - 901.0*e8/184320.0,
- 13.0*e4/768.0 + 17.0*e6/5120.0 - 311.0*e8/737280.0,
- 61.0*e6/15360.0 + 899.0*e8/430080.0,
- 49561.0*e8/41287680.0
- };
- double s = e * Math.sin(coord.lat());
- double l = Math.log(Math.tan(Math.PI/4.0 + coord.lat()/2.0) *
- Math.pow((1.0 - s) / (1.0 + s), e/2.0));
- double phi = Math.asin(Math.sin(coord.lon() - lg0) /
- ((Math.exp(l) + Math.exp(-l)) / 2.0));
- double ls = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
- double lambda = Math.atan(((Math.exp(l) - Math.exp(-l)) / 2.0) /
- Math.cos(coord.lon() - lg0));
-
- double north = C[0] * lambda;
- double east = C[0] * ls;
- for(int k = 1; k < 5; k++) {
- double r = 2.0 * k * lambda;
- double m = 2.0 * k * ls;
- double em = Math.exp(m);
- double en = Math.exp(-m);
- double sr = Math.sin(r)/2.0 * (em + en);
- double sm = Math.cos(r)/2.0 * (em - en);
- north += C[k] * sr;
- east += C[k] * sm;
- }
- east *= n;
- east += Xs;
- north *= n;
- north += Ys;
- return new EastNorth(east, north);
- }
-
- public LatLon eastNorth2latlon(EastNorth p) {
- if (currentGeodesic < 3) {
- MTProjection(p.east(), p.north(), zone, isNorth);
- LatLon geo = Geographic(p, Ellipsoid.hayford.a, Ellipsoid.hayford.e, 0.0 /* z */);
-
- // reference ellipsoid geographic => reference ellipsoid cartesian
- double N = Ellipsoid.hayford.a / (Math.sqrt(1.0 - Ellipsoid.hayford.e2 * Math.sin(geo.lat()) * Math.sin(geo.lat())));
- double X = (N /*+ h*/) * Math.cos(geo.lat()) * Math.cos(geo.lon());
- double Y = (N /*+ h*/) * Math.cos(geo.lat()) * Math.sin(geo.lon());
- double Z = (N * (1.0-Ellipsoid.hayford.e2) /*+ h*/) * Math.sin(geo.lat());
- // translation
- double coord[] = sevenParametersTransformation(X, Y, Z);
- // WGS84 cartesian => WGS84 geographic
- LatLon wgs = cart2LatLon(coord[0], coord[1], coord[2], Ellipsoid.GRS80);
- return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
- } else {
- // UTM_40S_Reunion_RGR92 or UTM_22N_Guyane_RGFG95
- LatLon geo = Geographic(p, Ellipsoid.GRS80.a, Ellipsoid.GRS80.e, 0.0 /* z */);
- double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(geo.lat()) * Math.sin(geo.lat())));
- double X = (N /*+ h*/) * Math.cos(geo.lat()) * Math.cos(geo.lon());
- double Y = (N /*+ h*/) * Math.cos(geo.lat()) * Math.sin(geo.lon());
- double Z = (N * (1.0-Ellipsoid.GRS80.e2) /*+ h*/) * Math.sin(geo.lat());
- LatLon wgs = cart2LatLon(X, Y, Z, Ellipsoid.GRS80);
- return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
- }
- }
-
- /**
- * initializes new projection coordinates (in north hemisphere)
- *
- * @param east east from origin in meters
- * @param north north from origin in meters
- * @param zone zone number (from 1 to 60)
- * @param isNorth true in north hemisphere, false in south hemisphere
- */
- private void MTProjection(double east, double north, int zone, boolean isNorth) {
- Ys = isNorth ? 0.0 : 10000000.0;
- double r6d = Math.PI / 30.0;
- lg0 = r6d * (zone - 0.5) - Math.PI;
- }
-
- public double scaleFactor() {
- return 1/Math.PI/2;
- }
-
- /**
- * initalizes from projected coordinates (Mercator transverse projection)
- *
- * @param coord projected coordinates pair
- * @param e reference ellipsoid excentricity
- * @param a reference ellipsoid long axis
- * @param z altitude in meters
- */
- private LatLon Geographic(EastNorth coord, double a, double e, double z) {
- double n = 0.9996 * a;
- double e2 = e * e;
- double e4 = e2 * e2;
- double e6 = e4 * e2;
- double e8 = e4 * e4;
- double C[] = {
- 1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
- e2/8.0 + e4/48.0 + 7.0*e6/2048.0 + e8/61440.0,
- e4/768.0 + 3.0*e6/1280.0 + 559.0*e8/368640.0,
- 17.0*e6/30720.0 + 283.0*e8/430080.0,
- 4397.0*e8/41287680.0
- };
- double l = (coord.north() - Ys) / (n * C[0]);
- double ls = (coord.east() - Xs) / (n * C[0]);
- double l0 = l;
- double ls0 = ls;
- for(int k = 1; k < 5; k++) {
- double r = 2.0 * k * l0;
- double m = 2.0 * k * ls0;
- double em = Math.exp(m);
- double en = Math.exp(-m);
- double sr = Math.sin(r)/2.0 * (em + en);
- double sm = Math.cos(r)/2.0 * (em - en);
- l -= C[k] * sr;
- ls -= C[k] * sm;
- }
- double lon = lg0 + Math.atan(((Math.exp(ls) - Math.exp(-ls)) / 2.0) /
- Math.cos(l));
- double phi = Math.asin(Math.sin(l) /
- ((Math.exp(ls) + Math.exp(-ls)) / 2.0));
- l = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
- double lat = 2.0 * Math.atan(Math.exp(l)) - Math.PI / 2.0;
- double lt0;
- do {
- lt0 = lat;
- double s = e * Math.sin(lat);
- lat = 2.0 * Math.atan(Math.pow((1.0 + s) / (1.0 - s), e/2.0) *
- Math.exp(l)) - Math.PI / 2.0;
- }
- while(Math.abs(lat-lt0) >= epsilon);
- //h = z;
-
- return new LatLon(lat, lon);
- }
-
- /**
- * initializes from cartesian coordinates
- *
- * @param X 1st coordinate in meters
- * @param Y 2nd coordinate in meters
- * @param Z 3rd coordinate in meters
- * @param ell reference ellipsoid
- */
- private LatLon cart2LatLon(double X, double Y, double Z, Ellipsoid ell) {
- double[] XYZ = {X, Y, Z};
- LatLon coord = ell.cart2LatLon(XYZ, epsilon);
- return new LatLon(Math.toRadians(coord.lat()), Math.toRadians(coord.lon()));
- }
-
- /**
- * 7 parameters transformation
- * @param coord X, Y, Z in array
- * @return transformed X, Y, Z in array
- */
- private double[] sevenParametersTransformation(double Xa, double Ya, double Za){
- double Xb = tx + Xa*(1+scaleDiff) + Za*ry - Ya*rz;
- double Yb = ty + Ya*(1+scaleDiff) + Xa*rz - Za*rx;
- double Zb = tz + Za*(1+scaleDiff) + Ya*rx - Xa*ry;
- return new double[]{Xb, Yb, Zb};
- }
-
- /**
- * 7 parameters inverse transformation
- * @param coord X, Y, Z in array
- * @return transformed X, Y, Z in array
- */
- private double[] invSevenParametersTransformation(double Xa, double Ya, double Za){
- double Xb = (1-scaleDiff)*(-tx + Xa + ((-tz+Za)*(-ry) - (-ty+Ya)*(-rz)));
- double Yb = (1-scaleDiff)*(-ty + Ya + ((-tx+Xa)*(-rz) - (-tz+Za)*(-rx)));
- double Zb = (1-scaleDiff)*(-tz + Za + ((-ty+Ya)*(-rx) - (-tx+Xa)*(-ry)));
- return new double[]{Xb, Yb, Zb};
- }
-
+ @Override
public String getCacheDirectoryName() {
return this.toString();
}
- /**
- * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
- */
- public double getDefaultZoomInPPD() {
- // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
- return 10.0;
- }
-
+ @Override
public Bounds getWorldBoundsLatLon() {
return utmBounds[currentGeodesic];
}
- public String toCode() {
+ @Override
+ public Integer getEpsgCode() {
return utmEPSGs[currentGeodesic];
}
@@ -401,12 +116,4 @@
public int hashCode() {
return getClass().getName().hashCode()+currentGeodesic; // our only real variable
- }
-
- @Override public String toString() {
- return (tr("UTM France (DOM)"));
- }
-
- public int getCurrentGeodesic() {
- return currentGeodesic;
}
@@ -426,4 +133,5 @@
}
+ @Override
public Collection getPreferences(JPanel p) {
Object prefcb = p.getComponent(2);
@@ -431,5 +139,4 @@
return null;
currentGeodesic = ((JComboBox)prefcb).getSelectedIndex();
- refresh7ParametersTranslation();
return Collections.singleton(Integer.toString(currentGeodesic+1));
}
@@ -437,16 +144,22 @@
@Override
public String[] allCodes() {
- return utmEPSGs;
+ String[] res = new String[utmEPSGs.length];
+ for (int i=0; i getPreferencesFromCode(String code) {
for (int i=0; i < utmEPSGs.length; i++ )
- if (utmEPSGs[i].endsWith(code))
+ if (("EPSG:"+utmEPSGs[i]).equals(code))
return Collections.singleton(Integer.toString(i+1));
return null;
}
+ @Override
public void setPreferences(Collection args) {
- currentGeodesic = DEFAULT_GEODESIC;
+ int currentGeodesic = DEFAULT_GEODESIC;
if (args != null) {
try {
@@ -461,6 +174,5 @@
} catch(NumberFormatException e) {}
}
- refresh7ParametersTranslation();
+ updateParameters(currentGeodesic);
}
-
}
Index: trunk/src/org/openstreetmap/josm/data/projection/datum/AbstractDatum.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/datum/AbstractDatum.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/datum/AbstractDatum.java (revision 4285)
@@ -0,0 +1,32 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.datum;
+
+import org.openstreetmap.josm.data.projection.Ellipsoid;
+
+abstract public class AbstractDatum implements Datum {
+
+ protected String name;
+ protected String proj4Id;
+ protected Ellipsoid ellps;
+
+ public AbstractDatum(String name, String proj4Id, Ellipsoid ellps) {
+ this.name = name;
+ this.proj4Id = proj4Id;
+ this.ellps = ellps;
+ }
+
+ @Override
+ public String getName() {
+ return name;
+ }
+
+ @Override
+ public String getProj4Id() {
+ return proj4Id;
+ }
+
+ @Override
+ public Ellipsoid getEllipsoid() {
+ return ellps;
+ }
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/datum/CentricDatum.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/datum/CentricDatum.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/datum/CentricDatum.java (revision 4285)
@@ -0,0 +1,27 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.datum;
+
+import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.Ellipsoid;
+
+/**
+ * A datum with different ellipsoid than WGS84, but does not require
+ * shift, rotation or scaling.
+ */
+public class CentricDatum extends AbstractDatum {
+
+ public CentricDatum(String name, String proj4Id, Ellipsoid ellps) {
+ super(name, proj4Id, ellps);
+ }
+
+ @Override
+ public LatLon toWGS84(LatLon ll) {
+ return Ellipsoid.WGS84.cart2LatLon(ellps.latLon2Cart(ll));
+ }
+
+ @Override
+ public LatLon fromWGS84(LatLon ll) {
+ return this.ellps.cart2LatLon(Ellipsoid.WGS84.latLon2Cart(ll));
+ }
+
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/datum/Datum.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/datum/Datum.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/datum/Datum.java (revision 4285)
@@ -0,0 +1,41 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.datum;
+
+import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.Ellipsoid;
+
+/**
+ * Represents a geodetic datum.
+ *
+ * Basically it provides conversion functions from and to the WGS84 datum.
+ */
+public interface Datum {
+
+ /**
+ * @return a human readable name of this projection
+ */
+ String getName();
+
+ /**
+ * @return the Proj.4 identifier
+ * (as reported by cs2cs -ld)
+ * If no id exists, return null.
+ */
+ String getProj4Id();
+
+ /**
+ * @return the ellipsoid associated with this datum
+ */
+ Ellipsoid getEllipsoid();
+
+ /**
+ * Convert lat/lon from this datum to WGS84 datum.
+ */
+ LatLon toWGS84(LatLon ll);
+
+ /**
+ * Convert lat/lon from WGS84 to this datum.
+ */
+ LatLon fromWGS84(LatLon ll);
+
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/datum/GRS80Datum.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/datum/GRS80Datum.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/datum/GRS80Datum.java (revision 4285)
@@ -0,0 +1,30 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.datum;
+
+import static org.openstreetmap.josm.tools.I18n.tr;
+
+import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.Ellipsoid;
+
+/**
+ * This datum indicates, that GRS80 ellipsoid is used and no conversion
+ * is necessary to get from or to the WGS84 datum.
+ */
+public class GRS80Datum extends AbstractDatum {
+
+ public static GRS80Datum INSTANCE = new GRS80Datum();
+
+ private GRS80Datum() {
+ super(tr("GRS80"), null, Ellipsoid.GRS80);
+ }
+
+ @Override
+ public LatLon fromWGS84(LatLon ll) {
+ return ll;
+ }
+
+ @Override
+ public LatLon toWGS84(LatLon ll) {
+ return ll;
+ }
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/datum/SevenParameterDatum.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/datum/SevenParameterDatum.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/datum/SevenParameterDatum.java (revision 4285)
@@ -0,0 +1,63 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.datum;
+
+import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.Ellipsoid;
+
+/**
+ * Datum provides general conversion from one ellipsoid to another.
+ *
+ * Seven parameters can be specified:
+ * - 3D offset
+ * - general rotation
+ * - scale
+ *
+ * This method is described by EPSG as EPSG::9606.
+ */
+public class SevenParameterDatum extends AbstractDatum {
+
+ protected double dx, dy, dz, rx, ry, rz, s;
+
+ /**
+ *
+ * @param name name of the datum
+ * @param proj4Id Proj.4 identifier for this datum (or null)
+ * @param ellps the ellipsoid used
+ * @param dx x offset in meters
+ * @param dy y offset in meters
+ * @param dz z offset in meters
+ * @param rx rotational parameter in seconds of arc
+ * @param ry rotational parameter in seconds of arc
+ * @param rz rotational parameter in seconds of arc
+ * @param s scale change in parts per million
+ */
+ public SevenParameterDatum(String name, String proj4Id, Ellipsoid ellps, double dx, double dy, double dz, double rx, double ry, double rz, double s) {
+ super(name, proj4Id, ellps);
+ this.dx = dx;
+ this.dy = dy;
+ this.dz = dz;
+ this.rx = Math.toRadians(rx / 3600);
+ this.ry = Math.toRadians(ry / 3600);
+ this.rz = Math.toRadians(rz / 3600);
+ this.s = s / 1e6;
+ }
+
+ @Override
+ public LatLon toWGS84(LatLon ll) {
+ double[] xyz = ellps.latLon2Cart(ll);
+ double x = dx + xyz[0]*(1+s) + xyz[2]*ry - xyz[1]*rz;
+ double y = dy + xyz[1]*(1+s) + xyz[0]*rz - xyz[2]*rx;
+ double z = dz + xyz[2]*(1+s) + xyz[1]*rx - xyz[0]*ry;
+ return Ellipsoid.WGS84.cart2LatLon(new double[] { x, y, z });
+ }
+
+ @Override
+ public LatLon fromWGS84(LatLon ll) {
+ double[] xyz = Ellipsoid.WGS84.latLon2Cart(ll);
+ double x = (1-s)*(-dx + xyz[0] + ((-dz+xyz[2])*(-ry) - (-dy+xyz[1])*(-rz)));
+ double y = (1-s)*(-dy + xyz[1] + ((-dx+xyz[0])*(-rz) - (-dz+xyz[2])*(-rx)));
+ double z = (1-s)*(-dz + xyz[2] + ((-dy+xyz[1])*(-rx) - (-dx+xyz[0])*(-ry)));
+ return this.ellps.cart2LatLon(new double[] { x, y, z });
+ }
+
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/datum/ThreeParameterDatum.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/datum/ThreeParameterDatum.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/datum/ThreeParameterDatum.java (revision 4285)
@@ -0,0 +1,39 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.datum;
+
+import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.Ellipsoid;
+
+/**
+ * Datum provides 3 dimensional offset and ellipsoid conversion.
+ */
+public class ThreeParameterDatum extends AbstractDatum {
+
+ protected double dx, dy, dz;
+
+ public ThreeParameterDatum(String name, String proj4Id, Ellipsoid ellps, double dx, double dy, double dz) {
+ super(name, proj4Id, ellps);
+ this.dx = dx;
+ this.dy = dy;
+ this.dz = dz;
+ }
+
+ @Override
+ public LatLon toWGS84(LatLon ll) {
+ double[] xyz = ellps.latLon2Cart(ll);
+ xyz[0] += dx;
+ xyz[1] += dy;
+ xyz[2] += dz;
+ return Ellipsoid.WGS84.cart2LatLon(xyz);
+ }
+
+ @Override
+ public LatLon fromWGS84(LatLon ll) {
+ double[] xyz = Ellipsoid.WGS84.latLon2Cart(ll);
+ xyz[0] -= dx;
+ xyz[1] -= dy;
+ xyz[2] -= dz;
+ return this.ellps.cart2LatLon(xyz);
+ }
+
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/datum/WGS84Datum.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/datum/WGS84Datum.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/datum/WGS84Datum.java (revision 4285)
@@ -0,0 +1,29 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.datum;
+
+import static org.openstreetmap.josm.tools.I18n.tr;
+
+import org.openstreetmap.josm.data.coor.LatLon;
+import org.openstreetmap.josm.data.projection.Ellipsoid;
+
+/**
+ * WGS84 datum. Transformation from and to WGS84 datum is a no-op.
+ */
+public class WGS84Datum extends AbstractDatum {
+
+ public static WGS84Datum INSTANCE = new WGS84Datum();
+
+ private WGS84Datum() {
+ super(tr("WGS84"), "WGS84", Ellipsoid.WGS84);
+ }
+
+ @Override
+ public LatLon fromWGS84(LatLon ll) {
+ return ll;
+ }
+
+ @Override
+ public LatLon toWGS84(LatLon ll) {
+ return ll;
+ }
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/proj/LambertConformalConic.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/proj/LambertConformalConic.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/proj/LambertConformalConic.java (revision 4285)
@@ -0,0 +1,151 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.proj;
+
+import static java.lang.Math.*;
+
+import static org.openstreetmap.josm.tools.I18n.tr;
+
+import org.openstreetmap.josm.data.projection.Ellipsoid;
+
+/**
+ * Implementation of the Lambert Conformal Conic projection.
+ *
+ * @author Pieren
+ */
+public class LambertConformalConic implements Proj {
+
+ protected Ellipsoid ellps;
+ protected double e;
+
+ /**
+ * projection exponent
+ */
+ protected double n;
+ /**
+ * projection factor
+ */
+ protected double F;
+ /**
+ * radius of the parallel of latitude of the false origin (2SP) or at
+ * natural origin (1SP)
+ */
+ protected double r0;
+
+ /**
+ * precision in iterative schema
+ */
+ protected static final double epsilon = 1e-12;
+
+ /**
+ * Constructor.
+ * Call one of the updateParameters... methods for initialization.
+ */
+ public LambertConformalConic() {
+ }
+
+ /**
+ * Initialize for LCC with 2 standard parallels.
+ *
+ * @param ellps the ellipsoid
+ * @param lat_0 latitude of false origin (in degrees)
+ * @param lat_1 latitude of first standard parallel (in degrees)
+ * @param lat_2 latitude of second standard parallel (in degrees)
+ */
+ public void updateParameters2SP(Ellipsoid ellps, double lat_0, double lat_1, double lat_2) {
+ this.ellps = ellps;
+ this.e = ellps.e;
+
+ final double m1 = m(toRadians(lat_1));
+ final double m2 = m(toRadians(lat_2));
+
+ final double t1 = t(toRadians(lat_1));
+ final double t2 = t(toRadians(lat_2));
+ final double tf = t(toRadians(lat_0));
+
+ n = (log(m1) - log(m2)) / (log(t1) - log(t2));
+ F = m1 / (n * pow(t1, n));
+ r0 = F * pow(tf, n);
+ }
+
+ /**
+ * Initialize for LCC with 1 standard parallel.
+ *
+ * @param ellps the ellipsoid
+ * @param lat_0 latitude of natural origin (in degrees)
+ */
+ public void updateParameters1SP(Ellipsoid ellps, double lat_0) {
+ this.ellps = ellps;
+ this.e = ellps.e;
+ final double lat_0_rad = toRadians(lat_0);
+
+ final double m0 = m(lat_0_rad);
+ final double t0 = t(lat_0_rad);
+
+ n = sin(lat_0_rad);
+ F = m0 / (n * pow(t0, n));
+ r0 = F * pow(t0, n);
+ }
+
+ /**
+ * Initialize LCC by providing the projection parameters directly.
+ *
+ * @param ellps the ellipsoid
+ * @param n see field n
+ * @param F see field F
+ * @param r0 see field r0
+ */
+ public void updateParametersDirect(Ellipsoid ellps, double n, double F, double r0) {
+ this.ellps = ellps;
+ this.e = ellps.e;
+ this.n = n;
+ this.F = F;
+ this.r0 = r0;
+ }
+
+ /**
+ * auxiliary function t
+ */
+ protected double t(double lat_rad) {
+ return tan(PI/4 - lat_rad / 2.0)
+ / pow(( (1.0 - e * sin(lat_rad)) / (1.0 + e * sin(lat_rad))) , e/2);
+ }
+
+ /**
+ * auxiliary function m
+ */
+ protected double m(double lat_rad) {
+ return cos(lat_rad) / (sqrt(1 - e * e * pow(sin(lat_rad), 2)));
+ }
+
+ @Override
+ public String getName() {
+ return tr("Lambert Conformal Conic");
+ }
+
+ @Override
+ public String getProj4Id() {
+ return "lcc";
+ }
+
+ @Override
+ public double[] project(double phi, double lambda) {
+ double sinphi = sin(phi);
+ double L = (0.5*log((1+sinphi)/(1-sinphi))) - e/2*log((1+e*sinphi)/(1-e*sinphi));
+ double r = F*exp(-n*L);
+ double gamma = n*lambda;
+ double X = r*sin(gamma);
+ double Y = r0 - r*cos(gamma);
+ return new double[] { X, Y };
+ }
+
+ @Override
+ public double[] invproject(double east, double north) {
+ double r = sqrt(pow(east,2) + pow(north-r0, 2));
+ double gamma = atan(east / (r0-north));
+ double lambda = gamma/n;
+ double latIso = (-1/n) * log(abs(r/F));
+ double phi = ellps.latitude(latIso, e, epsilon);
+ return new double[] { phi, lambda };
+ }
+
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/proj/Mercator.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/proj/Mercator.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/proj/Mercator.java (revision 4285)
@@ -0,0 +1,33 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.proj;
+
+import static java.lang.Math.*;
+
+import static org.openstreetmap.josm.tools.I18n.tr;
+
+/**
+ * Mercator Projection.
+ */
+public class Mercator implements Proj {
+
+ @Override
+ public String getName() {
+ return tr("Mercator");
+ }
+
+ @Override
+ public String getProj4Id() {
+ return "merc";
+ }
+
+ @Override
+ public double[] project(double lat_rad, double lon_rad) {
+ return new double[] { lon_rad, log(tan(PI/4 + lat_rad/2)) };
+ }
+
+ @Override
+ public double[] invproject(double east, double north) {
+ return new double[] { atan(sinh(north)), east };
+ }
+
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/proj/Proj.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/proj/Proj.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/proj/Proj.java (revision 4285)
@@ -0,0 +1,54 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.proj;
+
+/**
+ * A projection (in the narrow sense).
+ *
+ * Converts lat/lon the east/north and the other way around.
+ *
+ * Datum conversion, false easting / northing, origin of longitude
+ * and general scale factor is already applied when the projection is invoked.
+ *
+ * Lat/lon is not in degrees, but in radians (unlike other parts of JOSM).
+ * Additional parameters in the constructor arguments are usually still in
+ * degrees. So to avoid confusion, you can follow the convention, that
+ * coordinates in radians are called lat_rad/lon_rad or phi/lambda.
+ *
+ * East/north values are not in meters, but in meters divided by the semi major
+ * axis of the ellipsoid (earth radius). (Usually this is what you get anyway,
+ * unless you multiply by 'a' somehow implicitly or explicitly.)
+ *
+ */
+public interface Proj {
+ /**
+ * A Human readable name of this projection.
+ */
+ String getName();
+
+ /**
+ * The Proj.4 identifier.
+ *
+ * (as reported by cs2cs -lp)
+ * If no id exists, return null.
+ */
+ String getProj4Id();
+
+ /**
+ * Convert lat/lon to east/north.
+ *
+ * @param lat_rad the latitude in radians
+ * @param lon_rad the longitude in radians
+ * @return array of length 2, containing east and north value in meters,
+ * divided by the semi major axis of the ellipsoid.
+ */
+ double[] project(double lat_rad, double lon_rad);
+
+ /**
+ * Convert east/north to lat/lon.
+ *
+ * @param east east value in meters, divided by the semi major axis of the ellipsoid
+ * @param north north value in meters, divided by the semi major axis of the ellipsoid
+ * @return array of length 2, containing lat and lon in radians.
+ */
+ double[] invproject(double east, double north);
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/proj/SwissObliqueMercator.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/proj/SwissObliqueMercator.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/proj/SwissObliqueMercator.java (revision 4285)
@@ -0,0 +1,97 @@
+//License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.proj;
+
+import static java.lang.Math.*;
+
+import static org.openstreetmap.josm.tools.I18n.tr;
+
+import org.openstreetmap.josm.data.projection.Ellipsoid;
+
+/**
+ * Projection for the SwissGrid CH1903 / L03, see http://de.wikipedia.org/wiki/Swiss_Grid.
+ *
+ * Calculations are based on formula from
+ * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf
+ *
+ * August 2010 update to this formula (rigorous formulas)
+ * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf
+ */
+public class SwissObliqueMercator implements Proj {
+
+ private final Ellipsoid ellps;
+ private double kR;
+ private double alpha;
+ private double b0;
+ private double K;
+
+ private static final double EPSILON = 1e-11;
+
+ public SwissObliqueMercator(Ellipsoid ellps, double lat_0) {
+ this.ellps = ellps;
+ updateParameters(lat_0);
+ }
+
+ public void updateParameters(double lat_0) {
+ double phi0 = toRadians(lat_0);
+ kR = sqrt(1 - ellps.e2) / (1 - (ellps.e2 * pow(sin(phi0), 2)));
+ alpha = sqrt(1 + (ellps.eb2 * pow(cos(phi0), 4)));
+ b0 = asin(sin(phi0) / alpha);
+ K = log(tan(PI / 4 + b0 / 2)) - alpha
+ * log(tan(PI / 4 + phi0 / 2)) + alpha * ellps.e / 2
+ * log((1 + ellps.e * sin(phi0)) / (1 - ellps.e * sin(phi0)));
+ }
+
+ @Override
+ public String getName() {
+ return tr("Swiss Oblique Mercator");
+ }
+
+ @Override
+ public String getProj4Id() {
+ return "somerc";
+ }
+
+ @Override
+ public double[] project(double phi, double lambda) {
+
+ double S = alpha * log(tan(PI / 4 + phi / 2)) - alpha * ellps.e / 2
+ * log((1 + ellps.e * sin(phi)) / (1 - ellps.e * sin(phi))) + K;
+ double b = 2 * (atan(exp(S)) - PI / 4);
+ double l = alpha * lambda;
+
+ double lb = atan2(sin(l), sin(b0) * tan(b) + cos(b0) * cos(l));
+ double bb = asin(cos(b0) * sin(b) - sin(b0) * cos(b) * cos(l));
+
+ double y = kR * lb;
+ double x = kR / 2 * log((1 + sin(bb)) / (1 - sin(bb)));
+
+ return new double[] { y, x };
+ }
+
+ @Override
+ public double[] invproject(double y, double x) {
+ double lb = y / kR;
+ double bb = 2 * (atan(exp(x / kR)) - PI / 4);
+
+ double b = asin(cos(b0) * sin(bb) + sin(b0) * cos(bb) * cos(lb));
+ double l = atan2(sin(lb), cos(b0) * cos(lb) - sin(b0) * tan(bb));
+
+ double lambda = l / alpha;
+ double phi = b;
+ double S = 0;
+
+ double prevPhi = -1000;
+ int iteration = 0;
+ // iteration to finds S and phi
+ while (abs(phi - prevPhi) > EPSILON) {
+ if (++iteration > 30)
+ throw new RuntimeException("Two many iterations");
+ prevPhi = phi;
+ S = 1 / alpha * (log(tan(PI / 4 + b / 2)) - K) + ellps.e
+ * log(tan(PI / 4 + asin(ellps.e * sin(phi)) / 2));
+ phi = 2 * atan(exp(S)) - PI / 2;
+ }
+ return new double[] { phi, lambda };
+ }
+
+}
Index: trunk/src/org/openstreetmap/josm/data/projection/proj/TransverseMercator.java
===================================================================
--- trunk/src/org/openstreetmap/josm/data/projection/proj/TransverseMercator.java (revision 4285)
+++ trunk/src/org/openstreetmap/josm/data/projection/proj/TransverseMercator.java (revision 4285)
@@ -0,0 +1,283 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.proj;
+
+import static java.lang.Math.*;
+
+import static org.openstreetmap.josm.tools.I18n.tr;
+
+import org.openstreetmap.josm.data.projection.Ellipsoid;
+
+/**
+ * Transverse Mercator projection.
+ *
+ * @author Dirk Stöcker
+ * code based on JavaScript from Chuck Taylor
+ *
+ */
+public class TransverseMercator implements Proj {
+
+ protected double a, b;
+
+ public TransverseMercator(Ellipsoid ellps) {
+ this.a = ellps.a;
+ this.b = ellps.b;
+ }
+
+ @Override
+ public String getName() {
+ return tr("Transverse Mercator");
+ }
+
+ @Override
+ public String getProj4Id() {
+ return "tmerc";
+ }
+
+ /**
+ * Converts a latitude/longitude pair to x and y coordinates in the
+ * Transverse Mercator projection. Note that Transverse Mercator is not
+ * the same as UTM; a scale factor is required to convert between them.
+ *
+ * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
+ * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
+ *
+ * @param phi Latitude of the point, in radians
+ * @param lambda Longitude of the point, in radians
+ * @return A 2-element array containing the x and y coordinates
+ * of the computed point
+ */
+ @Override
+ public double[] project(double phi, double lambda) {
+
+ /* Precalculate ep2 */
+ double ep2 = (pow(a, 2.0) - pow(b, 2.0)) / pow(b, 2.0);
+
+ /* Precalculate nu2 */
+ double nu2 = ep2 * pow(cos(phi), 2.0);
+
+ /* Precalculate N / a */
+ double N_a = a / (b * sqrt(1 + nu2));
+
+ /* Precalculate t */
+ double t = tan(phi);
+ double t2 = t * t;
+
+ /* Precalculate l */
+ double l = lambda;
+
+ /* Precalculate coefficients for l**n in the equations below
+ so a normal human being can read the expressions for easting
+ and northing
+ -- l**1 and l**2 have coefficients of 1.0 */
+ double l3coef = 1.0 - t2 + nu2;
+
+ double l4coef = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2);
+
+ double l5coef = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2
+ - 58.0 * t2 * nu2;
+
+ double l6coef = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2
+ - 330.0 * t2 * nu2;
+
+ double l7coef = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2);
+
+ double l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
+
+ return new double[] {
+ /* Calculate easting (x) */
+ N_a * cos(phi) * l
+ + (N_a / 6.0 * pow(cos(phi), 3.0) * l3coef * pow(l, 3.0))
+ + (N_a / 120.0 * pow(cos(phi), 5.0) * l5coef * pow(l, 5.0))
+ + (N_a / 5040.0 * pow(cos(phi), 7.0) * l7coef * pow(l, 7.0)),
+ /* Calculate northing (y) */
+ ArcLengthOfMeridian (phi) / a
+ + (t / 2.0 * N_a * pow(cos(phi), 2.0) * pow(l, 2.0))
+ + (t / 24.0 * N_a * pow(cos(phi), 4.0) * l4coef * pow(l, 4.0))
+ + (t / 720.0 * N_a * pow(cos(phi), 6.0) * l6coef * pow(l, 6.0))
+ + (t / 40320.0 * N_a * pow(cos(phi), 8.0) * l8coef * pow(l, 8.0)) };
+ }
+
+ /**
+ * Converts x and y coordinates in the Transverse Mercator projection to
+ * a latitude/longitude pair. Note that Transverse Mercator is not
+ * the same as UTM; a scale factor is required to convert between them.
+ *
+ * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
+ * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
+ *
+ * Remarks:
+ * The local variables Nf, nuf2, tf, and tf2 serve the same purpose as
+ * N, nu2, t, and t2 in MapLatLonToXY, but they are computed with respect
+ * to the footpoint latitude phif.
+ *
+ * x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
+ * to optimize computations.
+ *
+ * @param x The easting of the point, in meters, divided by the semi major axis of the ellipsoid
+ * @param y The northing of the point, in meters, divided by the semi major axis of the ellipsoid
+ * @return A 2-element containing the latitude and longitude
+ * in radians
+ */
+ @Override
+ public double[] invproject(double x, double y) {
+ /* Get the value of phif, the footpoint latitude. */
+ double phif = footpointLatitude(y);
+
+ /* Precalculate ep2 */
+ double ep2 = (a*a - b*b)
+ / (b*b);
+
+ /* Precalculate cos(phif) */
+ double cf = cos(phif);
+
+ /* Precalculate nuf2 */
+ double nuf2 = ep2 * pow(cf, 2.0);
+
+ /* Precalculate Nf / a and initialize Nfpow */
+ double Nf_a = a / (b * sqrt(1 + nuf2));
+ double Nfpow = Nf_a;
+
+ /* Precalculate tf */
+ double tf = tan(phif);
+ double tf2 = tf * tf;
+ double tf4 = tf2 * tf2;
+
+ /* Precalculate fractional coefficients for x**n in the equations
+ below to simplify the expressions for latitude and longitude. */
+ double x1frac = 1.0 / (Nfpow * cf);
+
+ Nfpow *= Nf_a; /* now equals Nf**2) */
+ double x2frac = tf / (2.0 * Nfpow);
+
+ Nfpow *= Nf_a; /* now equals Nf**3) */
+ double x3frac = 1.0 / (6.0 * Nfpow * cf);
+
+ Nfpow *= Nf_a; /* now equals Nf**4) */
+ double x4frac = tf / (24.0 * Nfpow);
+
+ Nfpow *= Nf_a; /* now equals Nf**5) */
+ double x5frac = 1.0 / (120.0 * Nfpow * cf);
+
+ Nfpow *= Nf_a; /* now equals Nf**6) */
+ double x6frac = tf / (720.0 * Nfpow);
+
+ Nfpow *= Nf_a; /* now equals Nf**7) */
+ double x7frac = 1.0 / (5040.0 * Nfpow * cf);
+
+ Nfpow *= Nf_a; /* now equals Nf**8) */
+ double x8frac = tf / (40320.0 * Nfpow);
+
+ /* Precalculate polynomial coefficients for x**n.
+ -- x**1 does not have a polynomial coefficient. */
+ double x2poly = -1.0 - nuf2;
+ double x3poly = -1.0 - 2 * tf2 - nuf2;
+ double x4poly = 5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 - 3.0 * (nuf2 *nuf2) - 9.0 * tf2 * (nuf2 * nuf2);
+ double x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2;
+ double x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 + 162.0 * tf2 * nuf2;
+ double x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2);
+ double x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2);
+
+ return new double[] {
+ /* Calculate latitude */
+ phif + x2frac * x2poly * (x * x)
+ + x4frac * x4poly * pow(x, 4.0)
+ + x6frac * x6poly * pow(x, 6.0)
+ + x8frac * x8poly * pow(x, 8.0),
+ /* Calculate longitude */
+ x1frac * x
+ + x3frac * x3poly * pow(x, 3.0)
+ + x5frac * x5poly * pow(x, 5.0)
+ + x7frac * x7poly * pow(x, 7.0) };
+ }
+
+ /**
+ * ArcLengthOfMeridian
+ *
+ * Computes the ellipsoidal distance from the equator to a point at a
+ * given latitude.
+ *
+ * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
+ * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
+ *
+ * @param phi Latitude of the point, in radians
+ * @return The ellipsoidal distance of the point from the equator
+ * (in meters, divided by the semi major axis of the ellipsoid)
+ */
+ private double ArcLengthOfMeridian(double phi) {
+ /* Precalculate n */
+ double n = (a - b) / (a + b);
+
+ /* Precalculate alpha */
+ double alpha = ((a + b) / 2.0)
+ * (1.0 + (pow(n, 2.0) / 4.0) + (pow(n, 4.0) / 64.0));
+
+ /* Precalculate beta */
+ double beta = (-3.0 * n / 2.0) + (9.0 * pow(n, 3.0) / 16.0)
+ + (-3.0 * pow(n, 5.0) / 32.0);
+
+ /* Precalculate gamma */
+ double gamma = (15.0 * pow(n, 2.0) / 16.0)
+ + (-15.0 * pow(n, 4.0) / 32.0);
+
+ /* Precalculate delta */
+ double delta = (-35.0 * pow(n, 3.0) / 48.0)
+ + (105.0 * pow(n, 5.0) / 256.0);
+
+ /* Precalculate epsilon */
+ double epsilon = (315.0 * pow(n, 4.0) / 512.0);
+
+ /* Now calculate the sum of the series and return */
+ return alpha
+ * (phi + (beta * sin(2.0 * phi))
+ + (gamma * sin(4.0 * phi))
+ + (delta * sin(6.0 * phi))
+ + (epsilon * sin(8.0 * phi)));
+ }
+
+ /**
+ * FootpointLatitude
+ *
+ * Computes the footpoint latitude for use in converting transverse
+ * Mercator coordinates to ellipsoidal coordinates.
+ *
+ * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
+ * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
+ *
+ * @param y northing coordinate, in meters, divided by the semi major axis of the ellipsoid
+ * @return The footpoint latitude, in radians
+ */
+ private double footpointLatitude(double y) {
+ /* Precalculate n (Eq. 10.18) */
+ double n = (a - b) / (a + b);
+
+ /* Precalculate alpha_ (Eq. 10.22) */
+ /* (Same as alpha in Eq. 10.17) */
+ double alpha_ = ((a + b) / 2.0)
+ * (1 + (pow(n, 2.0) / 4) + (pow(n, 4.0) / 64));
+
+ /* Precalculate y_ (Eq. 10.23) */
+ double y_ = y / alpha_ * a;
+
+ /* Precalculate beta_ (Eq. 10.22) */
+ double beta_ = (3.0 * n / 2.0) + (-27.0 * pow(n, 3.0) / 32.0)
+ + (269.0 * pow(n, 5.0) / 512.0);
+
+ /* Precalculate gamma_ (Eq. 10.22) */
+ double gamma_ = (21.0 * pow(n, 2.0) / 16.0)
+ + (-55.0 * pow(n, 4.0) / 32.0);
+
+ /* Precalculate delta_ (Eq. 10.22) */
+ double delta_ = (151.0 * pow(n, 3.0) / 96.0)
+ + (-417.0 * pow(n, 5.0) / 128.0);
+
+ /* Precalculate epsilon_ (Eq. 10.22) */
+ double epsilon_ = (1097.0 * pow(n, 4.0) / 512.0);
+
+ /* Now calculate the sum of the series (Eq. 10.21) */
+ return y_ + (beta_ * sin(2.0 * y_))
+ + (gamma_ * sin(4.0 * y_))
+ + (delta_ * sin(6.0 * y_))
+ + (epsilon_ * sin(8.0 * y_));
+ }
+
+}