Index: /trunk/src/org/openstreetmap/josm/data/ProjectionBounds.java =================================================================== --- /trunk/src/org/openstreetmap/josm/data/ProjectionBounds.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/data/ProjectionBounds.java (revision 9419) @@ -123,4 +123,14 @@ /** + * Check, if a point is within the bounds. + * @param en the point + * @return true, if en is within the bounds + */ + public boolean contains(EastNorth en) { + return minEast <= en.east() && en.east() <= maxEast && + minNorth <= en.north() && en.north() <= maxNorth; + } + + /** * Returns the min east/north. * @return the min east/north Index: /trunk/src/org/openstreetmap/josm/data/coor/LatLon.java =================================================================== --- /trunk/src/org/openstreetmap/josm/data/coor/LatLon.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/data/coor/LatLon.java (revision 9419) @@ -56,4 +56,10 @@ public static final LatLon ZERO = new LatLon(0, 0); + /** + * North and south pole. + */ + public static final LatLon NORTH_POLE = new LatLon(90, 0); + public static final LatLon SOUTH_POLE = new LatLon(-90, 0); + private static DecimalFormat cDmsMinuteFormatter = new DecimalFormat("00"); private static DecimalFormat cDmsSecondFormatter = new DecimalFormat("00.0"); Index: /trunk/src/org/openstreetmap/josm/data/projection/CustomProjection.java =================================================================== --- /trunk/src/org/openstreetmap/josm/data/projection/CustomProjection.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/data/projection/CustomProjection.java (revision 9419) @@ -14,4 +14,6 @@ import org.openstreetmap.josm.Main; import org.openstreetmap.josm.data.Bounds; +import org.openstreetmap.josm.data.ProjectionBounds; +import org.openstreetmap.josm.data.coor.EastNorth; import org.openstreetmap.josm.data.coor.LatLon; import org.openstreetmap.josm.data.projection.datum.CentricDatum; @@ -23,4 +25,5 @@ import org.openstreetmap.josm.data.projection.datum.ThreeParameterDatum; import org.openstreetmap.josm.data.projection.datum.WGS84Datum; +import org.openstreetmap.josm.data.projection.proj.IPolar; import org.openstreetmap.josm.data.projection.proj.Mercator; import org.openstreetmap.josm.data.projection.proj.Proj; @@ -95,4 +98,6 @@ /** Latitude of second standard parallel */ lat_2("lat_2", true), + /** Latitude of true scale */ + lat_ts("lat_ts", true), /** the exact proj.4 string will be preserved in the WKT representation */ wktext("wktext", false), // ignored @@ -493,4 +498,8 @@ projParams.lat2 = parseAngle(s, Param.lat_2.key); } + s = parameters.get(Param.lat_ts.key); + if (s != null) { + projParams.lat_ts = parseAngle(s, Param.lat_ts.key); + } proj.initialize(projParams); return proj; @@ -690,3 +699,36 @@ return ret; } + + @Override + public Bounds getLatLonBoundsBox(ProjectionBounds r) { + Bounds result = new Bounds(eastNorth2latlon(r.getMin())); + result.extend(eastNorth2latlon(r.getMax())); + final int N = 40; + double dEast = (r.maxEast - r.minEast) / N; + double dNorth = (r.maxNorth - r.minNorth) / N; + for (int i = 0; i <= N; i++) { + result.extend(eastNorth2latlon(new EastNorth(r.minEast + i * dEast, r.minNorth))); + result.extend(eastNorth2latlon(new EastNorth(r.minEast + i * dEast, r.maxNorth))); + result.extend(eastNorth2latlon(new EastNorth(r.minEast, r.minNorth + i * dNorth))); + result.extend(eastNorth2latlon(new EastNorth(r.maxEast, r.minNorth + i * dNorth))); + } + // if the box contains one of the poles, the above method did not get + // correct min/max latitude value + if (proj instanceof IPolar) {; + IPolar polarProj = (IPolar) proj; + if (polarProj.hasPole(false)) { + EastNorth enNorthPole = latlon2eastNorth(LatLon.NORTH_POLE); + if (r.contains(enNorthPole)) { + result.extend(LatLon.NORTH_POLE); + } + } + if (polarProj.hasPole(true)) { + EastNorth enSouthPole = latlon2eastNorth(LatLon.SOUTH_POLE); + if (r.contains(enSouthPole)) { + result.extend(LatLon.SOUTH_POLE); + } + } + } + return result; + } } Index: /trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java =================================================================== --- /trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java (revision 9419) @@ -114,15 +114,18 @@ /** - * first eccentricity + * first eccentricity: + * sqrt(a*a - b*b) / a */ public final double e; /** - * first eccentricity squared + * first eccentricity squared: + * (a*a - b*b) / (a*a) */ public final double e2; /** - * square of the second eccentricity + * square of the second eccentricity: + * (a*a - b*b) / (b*b) */ public final double eb2; @@ -331,5 +334,5 @@ return new LatLon(Math.toDegrees(lt), Math.toDegrees(lg)); } - + /** * convert ellipsoidal coordinates to cartesian coordinates Index: /trunk/src/org/openstreetmap/josm/data/projection/Projection.java =================================================================== --- /trunk/src/org/openstreetmap/josm/data/projection/Projection.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/data/projection/Projection.java (revision 9419) @@ -83,4 +83,15 @@ /** + * Find lat/lon-box containing all the area of a given rectangle in + * east/north space. + * + * This is an approximate method. Points outside of the world should be ignored. + * + * @param pb the rectangle in projected space + * @return minimum lat/lon box, that when projected, covers pb + */ + Bounds getLatLonBoundsBox(ProjectionBounds pb); + + /** * Get the number of meters per unit of this projection. This more * defines the scale of the map, than real conversion of unit to meters Index: /trunk/src/org/openstreetmap/josm/data/projection/Projections.java =================================================================== --- /trunk/src/org/openstreetmap/josm/data/projection/Projections.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/data/projection/Projections.java (revision 9419) @@ -29,4 +29,5 @@ import org.openstreetmap.josm.data.projection.datum.ThreeParameterDatum; import org.openstreetmap.josm.data.projection.datum.WGS84Datum; +import org.openstreetmap.josm.data.projection.proj.AlbersEqualArea; import org.openstreetmap.josm.data.projection.proj.ClassProjFactory; import org.openstreetmap.josm.data.projection.proj.DoubleStereographic; @@ -34,4 +35,5 @@ import org.openstreetmap.josm.data.projection.proj.LonLat; import org.openstreetmap.josm.data.projection.proj.Mercator; +import org.openstreetmap.josm.data.projection.proj.PolarStereographic; import org.openstreetmap.josm.data.projection.proj.Proj; import org.openstreetmap.josm.data.projection.proj.ProjFactory; @@ -82,10 +84,12 @@ static { - registerBaseProjection("lonlat", LonLat.class, "core"); + registerBaseProjection("aea", AlbersEqualArea.class, "core"); registerBaseProjection("josm:smerc", Mercator.class, "core"); registerBaseProjection("lcc", LambertConformalConic.class, "core"); + registerBaseProjection("lonlat", LonLat.class, "core"); registerBaseProjection("somerc", SwissObliqueMercator.class, "core"); + registerBaseProjection("stere", PolarStereographic.class, "core"); + registerBaseProjection("sterea", DoubleStereographic.class, "core"); registerBaseProjection("tmerc", TransverseMercator.class, "core"); - registerBaseProjection("sterea", DoubleStereographic.class, "core"); ellipsoids.put("airy", Ellipsoid.Airy); Index: /trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java =================================================================== --- /trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java (revision 9419) @@ -53,4 +53,12 @@ /** + * Ellipsoid excentricity, equals to sqrt({@link #excentricitySquared}). + * Value 0 means that the ellipsoid is spherical. + * + * @see #excentricitySquared + */ + protected double e; + + /** * The square of excentricity: e² = (a²-b²)/a² where * e is the excentricity, @@ -63,4 +71,5 @@ public void initialize(ProjParameters params) throws ProjectionConfigurationException { e2 = params.ellps.e2; + e = params.ellps.e; // Compute constants for the mlfn double t; @@ -134,3 +143,25 @@ } } + + /** + * Computes function f(s,c,e²) = c/sqrt(1 - s²×e²) needed for the true scale + * latitude (Snyder 14-15), where s and c are the sine and cosine of + * the true scale latitude, and e² is the {@linkplain #excentricitySquared + * eccentricity squared}. + */ + final double msfn(final double s, final double c) { + return c / Math.sqrt(1.0 - (s*s) * e2); + } + + /** + * Computes function (15-9) and (9-13) from Snyder. + * Equivalent to negative of function (7-7). + */ + final double tsfn(final double lat, double sinlat) { + sinlat *= e; + /* + * NOTE: change sign to get the equivalent of Snyder (7-7). + */ + return Math.tan(0.5 * (Math.PI/2 - lat)) / Math.pow((1 - sinlat) / (1 + sinlat), 0.5*e); + } } Index: /trunk/src/org/openstreetmap/josm/data/projection/proj/AlbersEqualArea.java =================================================================== --- /trunk/src/org/openstreetmap/josm/data/projection/proj/AlbersEqualArea.java (revision 9419) +++ /trunk/src/org/openstreetmap/josm/data/projection/proj/AlbersEqualArea.java (revision 9419) @@ -0,0 +1,219 @@ +// License: GPL. For details, see LICENSE file. +package org.openstreetmap.josm.data.projection.proj; + +import static org.openstreetmap.josm.tools.I18n.tr; + +import org.openstreetmap.josm.data.Bounds; + +import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; + +/** + * Albers Equal Area Projection (EPSG code 9822). This is a conic projection with parallels being + * unequally spaced arcs of concentric circles, more closely spaced at north and south edges of the + * map. Merideans are equally spaced radii of the same circles and intersect parallels at right + * angles. As the name implies, this projection minimizes distortion in areas. + *

+ * The {@code "standard_parallel_2"} parameter is optional and will be given the same value as + * {@code "standard_parallel_1"} if not set (creating a 1 standard parallel projection). + *

+ * NOTE: + * formulae used below are from a port, to Java, of the {@code proj4} + * package of the USGS survey. USGS work is acknowledged here. + *

+ * This class has been derived from the implementation of the Geotools project; + * git 8cbf52d, org.geotools.referencing.operation.projection.AlbersEqualArea + * at the time of migration. + *

+* References: + *

+ * + * @see Albers Equal-Area Conic Projection on MathWorld + * @see "Albers_Conic_Equal_Area" on RemoteSensing.org + * @see British Columbia Albers Standard Projection + * + * @author Gerald I. Evenden (for original code in Proj4) + * @author Rueben Schulz + */ +public class AlbersEqualArea extends AbstractProj { + + /** + * Maximum number of iterations for iterative computations. + */ + private static final int MAXIMUM_ITERATIONS = 15; + + /** + * Difference allowed in iterative computations. + */ + private static final double ITERATION_TOLERANCE = 1E-10; + + /** + * Maximum difference allowed when comparing real numbers. + */ + private static final double EPSILON = 1E-6; + + /** + * Constants used by the spherical and elliptical Albers projection. + */ + private double n, c, rho0; + + /** + * An error condition indicating iteration will not converge for the + * inverse ellipse. See Snyder (14-20) + */ + private double ec; + + /** + * Standards parallel 1 in radians, for {@link #getParameterValues} implementation. + */ + private double phi1; + + /** + * Standards parallel 2 in radians, for {@link #getParameterValues} implementation. + */ + private double phi2; + + @Override + public String getName() { + return tr("Albers Equal Area"); + } + + @Override + public String getProj4Id() { + return "aea"; + } + + @Override + public void initialize(ProjParameters params) throws ProjectionConfigurationException { + super.initialize(params); + if (params.lat0 == null) + throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); + if (params.lat1 == null) + throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_1")); + + double lat0 = Math.toRadians(params.lat0); + phi1 = Math.toRadians(params.lat1); + phi2 = params.lat2 == null ? phi1 : Math.toRadians(params.lat2); + + // Compute Constants + if (Math.abs(phi1 + phi2) < EPSILON) { + throw new ProjectionConfigurationException(tr("standard parallels are opposite")); + } + double sinphi = Math.sin(phi1); + double cosphi = Math.cos(phi1); + double n = sinphi; + boolean secant = (Math.abs(phi1 - phi2) >= EPSILON); + double m1 = msfn(sinphi, cosphi); + double q1 = qsfn(sinphi); + if (secant) { // secant cone + sinphi = Math.sin(phi2); + cosphi = Math.cos(phi2); + double m2 = msfn(sinphi, cosphi); + double q2 = qsfn(sinphi); + n = (m1 * m1 - m2 * m2) / (q2 - q1); + } + c = m1 * m1 + n * q1; + rho0 = Math.sqrt(c - n * qsfn(Math.sin(lat0))) /n; + ec = 1.0 - .5 * (1.0-e2) * + Math.log((1.0 - e) / (1.0 + e)) / e; + this.n = n; + } + + @Override + public double[] project(double y, double x) { + x = normalizeLon(x); + x *= n; + double rho = c - n * qsfn(Math.sin(y)); + if (rho < 0.0) { + if (rho > -EPSILON) { + rho = 0.0; + } else { + throw new RuntimeException(); + } + } + rho = Math.sqrt(rho) / n; + y = rho0 - rho * Math.cos(x); + x = rho * Math.sin(x); + return new double[] {x, y}; + } + + @Override + public double[] invproject(double x, double y) { + y = rho0 - y; + double rho = Math.hypot(x, y); + if (rho > EPSILON) { + if (n < 0.0) { + rho = -rho; + x = -x; + y = -y; + } + x = Math.atan2(x, y) / n; + y = rho * n; + y = (c - y*y) / n; + if (Math.abs(y) <= ec) { + y = phi1(y); + } else { + y = (y < 0.0) ? -Math.PI/2.0 : Math.PI/2.0; + } + } else { + x = 0.0; + y = n > 0.0 ? Math.PI/2.0 : - Math.PI/2.0; + } + return new double[] {y, x}; + } + + /** + * Iteratively solves equation (3-16) from Snyder. + * + * @param qs arcsin(q/2), used in the first step of iteration + * @return the latitude + */ + public double phi1(final double qs) { + final double tone_es = 1 - e2; + double phi = Math.asin(0.5 * qs); + if (e < EPSILON) { + return phi; + } + for (int i=0; i + * + * In this file, only the polar case is implemented. This corresponds to + * EPSG:9810 (Polar Stereographic Variant A) and EPSG:9829 (Polar Stereographic + * Variant B). + *

+ * + * It is required, that the latitude of natural origin has the value +/-90 degrees. + *

+ * + * This class has been derived from the implementation of the Geotools project; + * git 8cbf52d, org.geotools.referencing.operation.projection.PolarStereographic + * at the time of migration. + *

+ * + * References: + *

+ * + * @see Stereographic projection on MathWorld + * @see Polar_Stereographic + * @see Oblique_Stereographic + * @see Stereographic + * @see Some Random Stereographic Issues + * + * @see DoubleStereographic + * + * @author André Gosselin + * @author Martin Desruisseaux (PMO, IRD) + * @author Rueben Schulz + */ +public class PolarStereographic extends AbstractProj implements IPolar { + /** + * Maximum number of iterations for iterative computations. + */ + private static final int MAXIMUM_ITERATIONS = 15; + + /** + * Difference allowed in iterative computations. + */ + private static final double ITERATION_TOLERANCE = 1E-10; + + /** + * Maximum difference allowed when comparing real numbers. + */ + private static final double EPSILON = 1E-8; + + /** + * A constant used in the transformations. + * This is not equal to the {@link #scaleFactor}. + */ + private double k0; + + /** + * Latitude of true scale, in radians + */ + private double latitudeTrueScale; + + /** + * {@code true} if this projection is for the south pole, or {@code false} + * if it is for the north pole. + */ + boolean southPole; + + @Override + public String getName() { + return tr("Polar Stereographic"); + } + + @Override + public String getProj4Id() { + return "stere"; + } + + @Override + public void initialize(ProjParameters params) throws ProjectionConfigurationException { + super.initialize(params); + if (params.lat0 == null) + throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); + if (params.lat0 != 90.0 && params.lat0 != -90.0) + throw new ProjectionConfigurationException( + tr("Polar Stereographic: Parameter ''{0}'' must be 90 or -90.", "lat_0")); + if (params.lat_ts == null) { + latitudeTrueScale = (params.lat0 < 0) ? -Math.PI/2 : Math.PI/2; + } else { + latitudeTrueScale = Math.toRadians(params.lat_ts); + } + southPole = latitudeTrueScale < 0; + + /* + * Computes coefficients. + */ + double latitudeTrueScaleAbs = Math.abs(latitudeTrueScale); + if (Math.abs(latitudeTrueScaleAbs - Math.PI/2) >= EPSILON) { + final double t = Math.sin(latitudeTrueScaleAbs); + k0 = msfn(t, Math.cos(latitudeTrueScaleAbs)) / + tsfn(latitudeTrueScaleAbs, t); // Derives from (21-32 and 21-33) + } else { + // True scale at pole (part of (21-33)) + k0 = 2.0 / Math.sqrt(Math.pow(1+e, 1+e)* + Math.pow(1-e, 1-e)); + } + } + + @Override + public double[] project(double y, double x) { + final double sinlat = Math.sin(y); + final double coslon = Math.cos(x); + final double sinlon = Math.sin(x); + if (southPole) { + final double rho = k0 * tsfn(-y, -sinlat); + x = rho * sinlon; + y = rho * coslon; + } else { + final double rho = k0 * tsfn(y, sinlat); + x = rho * sinlon; + y = -rho * coslon; + } + return new double[] {x, y}; + } + + @Override + public double[] invproject(double x, double y) { + final double rho = Math.hypot(x, y); + if (southPole) { + y = -y; + } + /* + * Compute latitude using iterative technique. + */ + final double t = rho/k0; + final double halfe = e/2.0; + double phi0 = 0; + for (int i=MAXIMUM_ITERATIONS;;) { + final double esinphi = e * Math.sin(phi0); + final double phi = (Math.PI/2) - 2.0*Math.atan(t*Math.pow((1-esinphi)/(1+esinphi), halfe)); + if (Math.abs(phi-phi0) < ITERATION_TOLERANCE) { + x = (Math.abs(rho) < EPSILON) ? 0.0 : Math.atan2(x, -y); + y = (southPole) ? -phi : phi; + break; + } + phi0 = phi; + if (--i < 0) { + throw new RuntimeException("no convergence"); + } + } + return new double[] {y, x}; + } + + @Override + public Bounds getAlgorithmBounds() { + final double CUT = 60; + if (southPole) { + return new Bounds(-90, -180, CUT, 180, false); + } else { + return new Bounds(-CUT, -180, 90, 180, false); + } + } + + @Override + public boolean hasPole(boolean south) { + return south == southPole; + } +} Index: /trunk/src/org/openstreetmap/josm/data/projection/proj/ProjParameters.java =================================================================== --- /trunk/src/org/openstreetmap/josm/data/projection/proj/ProjParameters.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/data/projection/proj/ProjParameters.java (revision 9419) @@ -14,3 +14,4 @@ public Double lat1; public Double lat2; + public Double lat_ts; } Index: /trunk/src/org/openstreetmap/josm/gui/NavigatableComponent.java =================================================================== --- /trunk/src/org/openstreetmap/josm/gui/NavigatableComponent.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/gui/NavigatableComponent.java (revision 9419) @@ -311,4 +311,12 @@ } + public ProjectionBounds getProjectionBounds(Rectangle r) { + EastNorth p1 = getEastNorth(r.x, r.y); + EastNorth p2 = getEastNorth(r.x + r.width, r.y + r.height); + ProjectionBounds pb = new ProjectionBounds(p1); + pb.extend(p2); + return pb; + } + /** * @param r rectangle @@ -316,25 +324,5 @@ */ public Bounds getLatLonBounds(Rectangle r) { - // TODO Maybe this should be (optional) method of Projection implementation - EastNorth p1 = getEastNorth(r.x, r.y); - EastNorth p2 = getEastNorth(r.x + r.width, r.y + r.height); - - Bounds result = new Bounds(Main.getProjection().eastNorth2latlon(p1)); - - double eastMin = Math.min(p1.east(), p2.east()); - double eastMax = Math.max(p1.east(), p2.east()); - double northMin = Math.min(p1.north(), p2.north()); - double northMax = Math.max(p1.north(), p2.north()); - double deltaEast = (eastMax - eastMin) / 10; - double deltaNorth = (northMax - northMin) / 10; - - for (int i = 0; i < 10; i++) { - result.extend(Main.getProjection().eastNorth2latlon(new EastNorth(eastMin + i * deltaEast, northMin))); - result.extend(Main.getProjection().eastNorth2latlon(new EastNorth(eastMin + i * deltaEast, northMax))); - result.extend(Main.getProjection().eastNorth2latlon(new EastNorth(eastMin, northMin + i * deltaNorth))); - result.extend(Main.getProjection().eastNorth2latlon(new EastNorth(eastMax, northMin + i * deltaNorth))); - } - - return result; + return Main.getProjection().getLatLonBoundsBox(getProjectionBounds(r)); } Index: /trunk/src/org/openstreetmap/josm/gui/util/HighlightHelper.java =================================================================== --- /trunk/src/org/openstreetmap/josm/gui/util/HighlightHelper.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/gui/util/HighlightHelper.java (revision 9419) @@ -14,6 +14,6 @@ /** - * This class stores the set of highlited primitives and - * allows easy and fast change of highlighting + * This class stores the set of highlighted primitives and + * allows easy and fast change of highlighting. */ public class HighlightHelper { Index: /trunk/src/org/openstreetmap/josm/tools/Utils.java =================================================================== --- /trunk/src/org/openstreetmap/josm/tools/Utils.java (revision 9418) +++ /trunk/src/org/openstreetmap/josm/tools/Utils.java (revision 9419) @@ -92,5 +92,5 @@ public static final String URL_CHARS = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789-._~:/?#[]@!$&'()*+,;=%"; - private static char[] DEFAULT_STRIP = {'\u200B', '\uFEFF'}; + private static final char[] DEFAULT_STRIP = {'\u200B', '\uFEFF'}; /** @@ -1000,5 +1000,6 @@ /** - * An alternative to {@link String#trim()} to effectively remove all leading and trailing white characters, including Unicode ones. + * An alternative to {@link String#trim()} to effectively remove all leading + * and trailing white characters, including Unicode ones. * @param str The string to strip * @return str, without leading and trailing characters, according to @@ -1017,5 +1018,6 @@ /** - * An alternative to {@link String#trim()} to effectively remove all leading and trailing white characters, including Unicode ones. + * An alternative to {@link String#trim()} to effectively remove all leading + * and trailing white characters, including Unicode ones. * @param str The string to strip * @param skipChars additional characters to skip