Index: trunk/src/org/openstreetmap/josm/data/projection/CustomProjection.java =================================================================== --- trunk/src/org/openstreetmap/josm/data/projection/CustomProjection.java (revision 13597) +++ trunk/src/org/openstreetmap/josm/data/projection/CustomProjection.java (revision 13598) @@ -296,7 +296,5 @@ } s = parameters.get(Param.bounds.key); - if (s != null) { - this.bounds = parseBounds(s); - } + this.bounds = s != null ? parseBounds(s) : null; s = parameters.get(Param.wmssrs.key); if (s != null) { @@ -590,4 +588,8 @@ projParams.gamma = parseAngle(s, Param.gamma.key); } + s = parameters.get(Param.lon_0.key); + if (s != null) { + projParams.lon0 = parseAngle(s, Param.lon_0.key); + } s = parameters.get(Param.lon_1.key); if (s != null) { Index: trunk/src/org/openstreetmap/josm/data/projection/Projections.java =================================================================== --- trunk/src/org/openstreetmap/josm/data/projection/Projections.java (revision 13597) +++ trunk/src/org/openstreetmap/josm/data/projection/Projections.java (revision 13598) @@ -25,7 +25,9 @@ import org.openstreetmap.josm.data.projection.datum.WGS84Datum; import org.openstreetmap.josm.data.projection.proj.AlbersEqualArea; +import org.openstreetmap.josm.data.projection.proj.AzimuthalEquidistant; import org.openstreetmap.josm.data.projection.proj.CassiniSoldner; import org.openstreetmap.josm.data.projection.proj.ClassProjFactory; import org.openstreetmap.josm.data.projection.proj.DoubleStereographic; +import org.openstreetmap.josm.data.projection.proj.EquidistantCylindrical; import org.openstreetmap.josm.data.projection.proj.LambertAzimuthalEqualArea; import org.openstreetmap.josm.data.projection.proj.LambertConformalConic; @@ -55,6 +57,15 @@ */ public static class ProjectionDefinition { + /** + * EPSG code + */ public final String code; + /** + * Projection name + */ public final String name; + /** + * projection definition (EPSG format) + */ public final String definition; @@ -89,5 +100,7 @@ static { registerBaseProjection("aea", AlbersEqualArea.class, "core"); + registerBaseProjection("aeqd", AzimuthalEquidistant.class, "core"); registerBaseProjection("cass", CassiniSoldner.class, "core"); + registerBaseProjection("eqc", EquidistantCylindrical.class, "core"); registerBaseProjection("laea", LambertAzimuthalEqualArea.class, "core"); registerBaseProjection("lcc", LambertConformalConic.class, "core"); @@ -199,4 +212,13 @@ } + /** + * Plugins can register additional base projections. + * + * @param id The "official" PROJ.4 id. In case the projection is not supported + * by PROJ.4, use some prefix, e.g. josm:myproj or gdal:otherproj. + * @param projClass The base projection class. + * @param origin Multiple plugins may implement the same base projection. + * Provide plugin name or similar string, so it be differentiated. + */ public static void registerBaseProjection(String id, Class projClass, String origin) { registerBaseProjection(id, new ClassProjFactory(projClass), origin); @@ -293,21 +315,30 @@ List result = new ArrayList<>(); Pattern epsgPattern = Pattern.compile("<(\\d+)>(.*)<>"); - String line, lastline = ""; + StringBuilder sb = new StringBuilder(); + String line; while ((line = r.readLine()) != null) { line = line.trim(); - if (!line.startsWith("#") && !line.isEmpty()) { - if (!lastline.startsWith("#")) throw new AssertionError("EPSG file seems corrupted"); - String name = lastline.substring(1).trim(); - Matcher m = epsgPattern.matcher(line); - if (m.matches()) { - String code = "EPSG:" + m.group(1); - String definition = m.group(2).trim(); - result.add(new ProjectionDefinition(code, name, definition)); - } else { - Logging.warn("Failed to parse line from the EPSG projection definition: "+line); + if (!line.isEmpty()) { + if (!line.startsWith("#")) { + Matcher m = epsgPattern.matcher(line); + if (m.matches()) { + String code = "EPSG:" + m.group(1); + String definition = m.group(2).trim(); + result.add(new ProjectionDefinition(code, sb.toString(), definition)); + } else { + Logging.warn("Failed to parse line from the EPSG projection definition: "+line); + } + sb.setLength(0); + } else if (!line.startsWith("# area: ")) { + if (sb.length() == 0) { + sb.append(line.substring(1).trim()); + } else { + sb.append('(').append(line.substring(1).trim()).append(')'); + } } } - lastline = line; - } + } + if (result.isEmpty()) + throw new AssertionError("EPSG file seems corrupted"); return result; } Index: trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java =================================================================== --- trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java (revision 13597) +++ trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java (revision 13598) @@ -155,4 +155,17 @@ } + /** + * Tolerant asin that will just return the limits of its output range if the input is out of range + * @param v the value whose arc sine is to be returned. + * @return the arc sine of the argument. + */ + protected final double aasin(double v) { + double av = Math.abs(v); + if (av >= 1.) { + return (v < 0. ? -Math.PI / 2 : Math.PI / 2); + } + return Math.asin(v); + } + // Iteratively solve equation (7-9) from Snyder. final double cphi2(final double ts) { Index: trunk/src/org/openstreetmap/josm/data/projection/proj/AzimuthalEquidistant.java =================================================================== --- trunk/src/org/openstreetmap/josm/data/projection/proj/AzimuthalEquidistant.java (revision 13598) +++ trunk/src/org/openstreetmap/josm/data/projection/proj/AzimuthalEquidistant.java (revision 13598) @@ -0,0 +1,308 @@ +// 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; +import org.openstreetmap.josm.tools.JosmRuntimeException; + +/** + * Azimuthal Equidistant projection. + *

+ * This implementation does not include the Guam or Micronesia variants. + * + * @author Gerald Evenden (original PROJ.4 implementation in C) + * @author Ben Caradoc-Davies (Transient Software Limited) + * @see Map Projections: A Working Manual, Snyder (1987), pages 191-202 + * @see PROJ.4 notes on parameters + * @see PROJ.4 implemention in C + * @see Wikipedia + * @see Wolfram Alpha + * @since 13598 + */ +public class AzimuthalEquidistant extends AbstractProj { + + /** + * Less strict tolerance. + */ + public static final double EPS10 = 1.e-10; + + /** + * Stricter tolerance. + */ + public static final double TOL = 1.e-14; + + /** + * Half of π. + */ + public static final double HALF_PI = Math.PI / 2; + + /** + * The four possible modes or aspects of the projection. + */ + public enum Mode { + /** North pole */ + NORTH_POLAR, + /** South pole */ + SOUTH_POLAR, + /** Equator */ + EQUATORIAL, + /** Oblique */ + OBLIQUE; + } + + /** + * Length of semi-major axis, in metres. This is named 'a' or 'R' + * (Radius in spherical cases) in Snyder. + */ + protected double semiMajor; + + /** + * Length of semi-minor axis, in metres. This is named 'b' in Snyder. + */ + protected double semiMinor; + + /** + * Central longitude in radians. Default value is 0, the Greenwich meridian. + * This is called 'lambda0' in Snyder. + */ + protected double centralMeridian; + + /** + * Latitude of origin in radians. Default value is 0, the equator. + * This is called 'phi0' in Snyder. + */ + protected double latitudeOfOrigin; + + /** + * The mode or aspect of the projection. + */ + protected Mode mode; + + /** + * Geodesic calculator used for this projection. Not used and set to null for polar projections. + */ + //protected Geodesic geodesic; // See https://josm.openstreetmap.de/ticket/16129#comment:21 + + /** + * The sine of the central latitude of the projection. + */ + protected double sinph0; + + /** + * The cosine of the central latitude of the projection. + */ + protected double cosph0; + + /** + * Meridian distance from the equator to the pole. Not used and set to NaN for non-polar projections. + */ + protected double mp; + + @Override + public String getName() { + return tr("Azimuthal Equidistant"); + } + + @Override + public String getProj4Id() { + return "aeqd"; + } + + @Override + public void initialize(ProjParameters params) throws ProjectionConfigurationException { + super.initialize(params); + if (params.lon0 == null) + throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lon_0")); + if (params.lat0 == null) + throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); + centralMeridian = Math.toRadians(params.lon0); + latitudeOfOrigin = Math.toRadians(params.lat0); + semiMajor = params.ellps.a; + semiMinor = params.ellps.b; + if (Math.abs(latitudeOfOrigin - HALF_PI) < EPS10) { + mode = Mode.NORTH_POLAR; + mp = mlfn(HALF_PI, 1, 0); + sinph0 = 1; + cosph0 = 0; + } else if (Math.abs(latitudeOfOrigin + HALF_PI) < EPS10) { + mode = Mode.SOUTH_POLAR; + mp = mlfn(-HALF_PI, -1, 0); + sinph0 = -1; + cosph0 = 0; + } else if (Math.abs(latitudeOfOrigin) < EPS10) { + mode = Mode.EQUATORIAL; + mp = Double.NaN; + sinph0 = 0; + cosph0 = 1; + //geodesic = new Geodesic(semiMajor, (semiMajor - semiMinor) / semiMajor); + throw new ProjectionConfigurationException("Equatorial AzimuthalEquidistant not yet supported"); + } else { + mode = Mode.OBLIQUE; + mp = Double.NaN; + sinph0 = Math.sin(latitudeOfOrigin); + cosph0 = Math.cos(latitudeOfOrigin); + //geodesic = new Geodesic(semiMajor, (semiMajor - semiMinor) / semiMajor); + throw new ProjectionConfigurationException("Oblique AzimuthalEquidistant not yet supported"); + } + } + + @Override + public Bounds getAlgorithmBounds() { + return new Bounds(-89, -180, 89, 180, false); + } + + @Override + public double[] project(double latRad, double lonRad) { + return spherical ? projectSpherical(latRad, lonRad) : projectEllipsoidal(latRad, lonRad); + } + + @Override + public double[] invproject(double east, double north) { + return spherical ? invprojectSpherical(east, north) : invprojectEllipsoidal(east, north); + } + + double[] projectSpherical(double latRad, double lonRad) { + double x = 0; + double y = 0; + double sinphi = Math.sin(latRad); + double cosphi = Math.cos(latRad); + double coslam = Math.cos(lonRad); + switch (mode) { + case EQUATORIAL: + case OBLIQUE: + if (mode == Mode.EQUATORIAL) { + y = cosphi * coslam; + } else { // Oblique + y = sinph0 * sinphi + cosph0 * cosphi * coslam; + } + if (Math.abs(Math.abs(y) - 1) < TOL) { + if (y < 0) { + throw new JosmRuntimeException("TOLERANCE_ERROR"); + } else { + x = 0; + y = 0; + } + } else { + y = Math.acos(y); + y /= Math.sin(y); + x = y * cosphi * Math.sin(lonRad); + y *= (mode == Mode.EQUATORIAL) ? sinphi + : (cosph0 * sinphi - sinph0 * cosphi * coslam); + } + break; + case NORTH_POLAR: + latRad = -latRad; + coslam = -coslam; + case SOUTH_POLAR: + if (Math.abs(latRad - HALF_PI) < EPS10) { + throw new JosmRuntimeException("TOLERANCE_ERROR"); + } + y = HALF_PI + latRad; + x = y * Math.sin(lonRad); + y *= coslam; + break; + } + return new double[] {x, y}; + } + + double[] invprojectSpherical(double east, double north) { + double x = east; + double y = north; + double lambda = 0; + double phi = 0; + double c_rh = Math.hypot(x, y); + if (c_rh > Math.PI) { + if (c_rh - EPS10 > Math.PI) { + throw new JosmRuntimeException("TOLERANCE_ERROR"); + } + c_rh = Math.PI; + } else if (c_rh < EPS10) { + phi = latitudeOfOrigin; + lambda = 0.; + } else { + if (mode == Mode.OBLIQUE || mode == Mode.EQUATORIAL) { + double sinc = Math.sin(c_rh); + double cosc = Math.cos(c_rh); + if (mode == Mode.EQUATORIAL) { + phi = aasin(y * sinc / c_rh); + x *= sinc; + y = cosc * c_rh; + } else { // Oblique + phi = aasin(cosc * sinph0 + y * sinc * cosph0 / c_rh); + y = (cosc - sinph0 * Math.sin(phi)) * c_rh; + x *= sinc * cosph0; + } + lambda = (y == 0) ? 0 : Math.atan2(x, y); + } else if (mode == Mode.NORTH_POLAR) { + phi = HALF_PI - c_rh; + lambda = Math.atan2(x, -y); + } else { // South Polar + phi = c_rh - HALF_PI; + lambda = Math.atan2(x, y); + } + } + return new double[] {phi, lambda}; + } + + double[] projectEllipsoidal(double latRad, double lonRad) { + double x = 0; + double y = 0; + double coslam = Math.cos(lonRad); + double cosphi = Math.cos(latRad); + double sinphi = Math.sin(latRad); + switch (mode) { + case NORTH_POLAR: + coslam = -coslam; + case SOUTH_POLAR: + double rho = Math.abs(mp - mlfn(latRad, sinphi, cosphi)); + x = rho * Math.sin(lonRad); + y = rho * coslam; + break; + case EQUATORIAL: + case OBLIQUE: + if (Math.abs(lonRad) < EPS10 && Math.abs(latRad - latitudeOfOrigin) < EPS10) { + x = 0; + y = 0; + break; + } + /*GeodesicData g = geodesic.Inverse(Math.toDegrees(latitudeOfOrigin), + Math.toDegrees(centralMeridian), Math.toDegrees(latRad), + Math.toDegrees(lonRad + centralMeridian)); + double azi1 = Math.toRadians(g.azi1); + x = g.s12 * Math.sin(azi1) / semiMajor; + y = g.s12 * Math.cos(azi1) / semiMajor;*/ + break; + } + return new double[] {x, y}; + } + + double[] invprojectEllipsoidal(double east, double north) { + double x = east; + double y = north; + double lambda = 0; + double phi = 0; + double c = Math.hypot(x, y); + if (c < EPS10) { + phi = latitudeOfOrigin; + lambda = 0; + } else { + if (mode == Mode.OBLIQUE || mode == Mode.EQUATORIAL) { + /*double x2 = x * semiMajor; + double y2 = y * semiMajor; + double azi1 = Math.atan2(x2, y2); + double s12 = Math.sqrt(x2 * x2 + y2 * y2); + GeodesicData g = geodesic.Direct(Math.toDegrees(latitudeOfOrigin), + Math.toDegrees(centralMeridian), Math.toDegrees(azi1), s12); + phi = Math.toRadians(g.lat2); + lambda = Math.toRadians(g.lon2);*/ + lambda -= centralMeridian; + } else { // Polar + phi = invMlfn((mode == Mode.NORTH_POLAR) ? (mp - c) : (mp + c)); + lambda = Math.atan2(x, (mode == Mode.NORTH_POLAR) ? -y : y); + } + } + return new double[] {phi, lambda}; + } +} Index: trunk/src/org/openstreetmap/josm/data/projection/proj/EquidistantCylindrical.java =================================================================== --- trunk/src/org/openstreetmap/josm/data/projection/proj/EquidistantCylindrical.java (revision 13598) +++ trunk/src/org/openstreetmap/josm/data/projection/proj/EquidistantCylindrical.java (revision 13598) @@ -0,0 +1,79 @@ +// 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; + +/** + * Equidistant cylindrical projection (EPSG code 9823). + * In the particular case where the {@code standard_parallel_1} is 0°, this projection is also called Plate Carree or Equirectangular. + * This is used in, for example, WGS84 / Plate Carree (EPSG:32662). + *

+ * References: + *

+ * + * @author John Grange + * @author Martin Desruisseaux + * + * @see Cylindrical Equidistant projection on MathWorld + * @see "Equirectangular" on RemoteSensing.org + * @since 13598 + */ +public class EquidistantCylindrical extends AbstractProj { + + /** + * Cosinus of the {@code "standard_parallel_1"} parameter. + */ + private double cosStandardParallel; + + /** + * Standard parallel parameter. + * Set to 0° for the {@code PlateCarree} case. + */ + private double standardParallel; + + @Override + public String getName() { + return tr("Equidistant Cylindrical (Plate Caree)"); + } + + @Override + public String getProj4Id() { + return "eqc"; + } + + @Override + public void initialize(ProjParameters params) throws ProjectionConfigurationException { + super.initialize(params); + if (params.lat1 != null) { + standardParallel = Math.abs(params.lat1); + cosStandardParallel = Math.cos(standardParallel); + } else { + // standard parallel is the equator (Plate Carree or Equirectangular) + standardParallel = 0; + cosStandardParallel = 1.0; + } + } + + @Override + public double[] project(double latRad, double lonRad) { + return new double[] {lonRad * cosStandardParallel, latRad}; + } + + @Override + public double[] invproject(double east, double north) { + return new double[] {north, east / cosStandardParallel}; + } + + @Override + public Bounds getAlgorithmBounds() { + return new Bounds(-89, -180, 89, 180, false); + } +} Index: trunk/src/org/openstreetmap/josm/data/projection/proj/ProjParameters.java =================================================================== --- trunk/src/org/openstreetmap/josm/data/projection/proj/ProjParameters.java (revision 13597) +++ trunk/src/org/openstreetmap/josm/data/projection/proj/ProjParameters.java (revision 13598) @@ -2,26 +2,43 @@ package org.openstreetmap.josm.data.projection.proj; +import org.openstreetmap.josm.data.projection.CustomProjection.Param; import org.openstreetmap.josm.data.projection.Ellipsoid; /** * Parameters to initialize a Proj object. + * @since 5066 */ public class ProjParameters { + /** {@code +ellps} */ public Ellipsoid ellps; + /** {@link Param#lat_0} */ public Double lat0; + /** {@link Param#lat_1} */ public Double lat1; + /** {@link Param#lat_2} */ public Double lat2; // Polar Stereographic and Mercator + /** {@link Param#lat_ts} */ public Double lat_ts; + // Azimuthal Equidistant + /** {@link Param#lon_0} */ + public Double lon0; + // Oblique Mercator + /** {@link Param#lonc} */ public Double lonc; + /** {@link Param#alpha} */ public Double alpha; + /** {@link Param#gamma} */ public Double gamma; + /** {@link Param#no_off} */ public Boolean no_off; + /** {@link Param#lon_1} */ public Double lon1; + /** {@link Param#lon_2} */ public Double lon2; }