Index: /trunk/CONTRIBUTION =================================================================== --- /trunk/CONTRIBUTION (revision 13597) +++ /trunk/CONTRIBUTION (revision 13598) @@ -93,4 +93,6 @@ Swiss CHENYX06 NTV2 grid: Source: Swiss Federal Office of Topography +ESRI projection definitions: Environmental Systems Research Institute + ------------------------------------ ICONS ------------------------------------ Index: /trunk/scripts/BuildProjectionDefinitions.java =================================================================== --- /trunk/scripts/BuildProjectionDefinitions.java (revision 13597) +++ /trunk/scripts/BuildProjectionDefinitions.java (revision 13598) @@ -7,9 +7,13 @@ import java.io.OutputStreamWriter; import java.nio.charset.StandardCharsets; -import java.util.HashMap; +import java.util.Arrays; import java.util.LinkedHashMap; +import java.util.List; +import java.util.Locale; import java.util.Map; +import java.util.TreeMap; import org.openstreetmap.josm.data.projection.CustomProjection; +import org.openstreetmap.josm.data.projection.CustomProjection.Param; import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; import org.openstreetmap.josm.data.projection.Projections; @@ -23,7 +27,8 @@ public class BuildProjectionDefinitions { - private static final String JOSM_EPSG_FILE = "data_nodist/projection/josm-epsg"; - private static final String PROJ4_EPSG_FILE = "data_nodist/projection/epsg"; - private static final String PROJ4_ESRI_FILE = "data_nodist/projection/esri"; + private static final String PROJ_DIR = "data_nodist/projection"; + private static final String JOSM_EPSG_FILE = "josm-epsg"; + private static final String PROJ4_EPSG_FILE = "epsg"; + private static final String PROJ4_ESRI_FILE = "esri"; private static final String OUTPUT_EPSG_FILE = "data/projection/custom-epsg"; @@ -40,5 +45,6 @@ private static int noGeocent = 0; private static int noBaseProjection = 0; - private static final Map baseProjectionMap = new HashMap<>(); + private static int noEllipsoid = 0; + private static int noNadgrid = 0; private static int noDatumgrid = 0; private static int noJosm = 0; @@ -47,4 +53,12 @@ private static int noOmercNoBounds = 0; private static int noEquatorStereo = 0; + + private static final Map baseProjectionMap = new TreeMap<>(); + private static final Map ellipsoidMap = new TreeMap<>(); + private static final Map nadgridMap = new TreeMap<>(); + private static final Map datumgridMap = new TreeMap<>(); + + private static List knownGeoidgrids; + private static List knownNadgrids; /** @@ -57,6 +71,12 @@ } + static List initList(String baseDir, String ext) { + return Arrays.asList(new File(baseDir + File.separator + PROJ_DIR) + .list((dir, name) -> !name.contains(".") || name.toLowerCase(Locale.ENGLISH).endsWith(ext))); + } + static void initMap(String baseDir, String file, Map map) throws IOException { - for (ProjectionDefinition pd : Projections.loadProjectionDefinitions(baseDir + File.separator + file)) { + for (ProjectionDefinition pd : Projections.loadProjectionDefinitions( + baseDir + File.separator + PROJ_DIR + File.separator + file)) { map.put(pd.code, pd); } @@ -67,4 +87,7 @@ initMap(baseDir, PROJ4_EPSG_FILE, epsgProj4); initMap(baseDir, PROJ4_ESRI_FILE, esriProj4); + + knownGeoidgrids = initList(baseDir, ".gtx"); + knownNadgrids = initList(baseDir, ".gsb"); try (FileOutputStream output = new FileOutputStream(baseDir + File.separator + OUTPUT_EPSG_FILE); @@ -104,10 +127,25 @@ System.out.println("some entries from proj.4 have not been included:"); System.out.println(String.format(" * already in the maintained JOSM list: %d entries", noInJosm)); - System.out.println(String.format(" * ESRI already in the standard EPSG list: %d entries", noInProj4)); + if (noInProj4 > 0) { + System.out.println(String.format(" * ESRI already in the standard EPSG list: %d entries", noInProj4)); + } System.out.println(String.format(" * deprecated: %d entries", noDeprecated)); System.out.println(String.format(" * using +proj=geocent, which is 3D (X,Y,Z) and not useful in JOSM: %d entries", noGeocent)); - System.out.println(String.format(" * unsupported base projection: %d entries", noBaseProjection)); - System.out.println(" in particular: " + baseProjectionMap); - System.out.println(String.format(" * requires data file for datum conversion: %d entries", noDatumgrid)); + if (noEllipsoid > 0) { + System.out.println(String.format(" * unsupported ellipsoids: %d entries", noEllipsoid)); + System.out.println(" in particular: " + ellipsoidMap); + } + if (noBaseProjection > 0) { + System.out.println(String.format(" * unsupported base projection: %d entries", noBaseProjection)); + System.out.println(" in particular: " + baseProjectionMap); + } + if (noDatumgrid > 0) { + System.out.println(String.format(" * requires data file for vertical datum conversion: %d entries", noDatumgrid)); + System.out.println(" in particular: " + datumgridMap); + } + if (noNadgrid > 0) { + System.out.println(String.format(" * requires data file for datum conversion: %d entries", noNadgrid)); + System.out.println(" in particular: " + nadgridMap); + } if (noOmercNoBounds > 0) { System.out.println(String.format(" * projection is Oblique Mercator (requires bounds), but no bounds specified: %d entries", noOmercNoBounds)); @@ -147,7 +185,8 @@ } - // exclude deprecated projections + // exclude deprecated/discontinued projections // EPSG:4296 is also deprecated, but this is not mentioned in the name - if (pd.name.contains("deprecated") || pd.code.equals("EPSG:4296")) { + String lowName = pd.name.toLowerCase(Locale.ENGLISH); + if (lowName.contains("deprecated") || lowName.contains("discontinued") || pd.code.equals("EPSG:4296")) { result = false; noDeprecated++; @@ -176,8 +215,18 @@ } - // requires datum conversion database - if (parameters.containsKey("geoidgrids")) { + // requires vertical datum conversion database (.gtx) + String geoidgrids = parameters.get("geoidgrids"); + if (geoidgrids != null && !"@null".equals(geoidgrids) && !knownGeoidgrids.contains(geoidgrids)) { result = false; noDatumgrid++; + incMap(datumgridMap, geoidgrids); + } + + // requires datum conversion database (.gsb) + String nadgrids = parameters.get("nadgrids"); + if (nadgrids != null && !"@null".equals(nadgrids) && !knownNadgrids.contains(nadgrids)) { + result = false; + noNadgrid++; + incMap(nadgridMap, nadgrids); } @@ -188,22 +237,65 @@ noBaseProjection++; if (!"geocent".equals(proj)) { - if (!baseProjectionMap.containsKey(proj)) { - baseProjectionMap.put(proj, 0); + incMap(baseProjectionMap, proj); + } + } + + // exclude entries where we don't support the base ellipsoid + String ellps = parameters.get("ellps"); + if (result && ellps != null && Projections.getEllipsoid(ellps) == null) { + result = false; + noEllipsoid++; + incMap(ellipsoidMap, ellps); + } + + if (result && "omerc".equals(proj) && !parameters.containsKey(CustomProjection.Param.bounds.key)) { + result = false; + noOmercNoBounds++; + } + + final double EPS10 = 1.e-10; + + String lat0 = parameters.get("lat_0"); + if (lat0 != null) { + try { + final double latitudeOfOrigin = Math.toRadians(CustomProjection.parseAngle(lat0, Param.lat_0.key)); + // TODO: implement equatorial stereographic, see https://josm.openstreetmap.de/ticket/15970 + if (result && "stere".equals(proj) && Math.abs(latitudeOfOrigin) < EPS10) { + result = false; + noEquatorStereo++; } - baseProjectionMap.put(proj, baseProjectionMap.get(proj)+1); - } - } - - if (result && "omerc".equals(proj) && !parameters.containsKey(CustomProjection.Param.bounds.key)) { - result = false; - noOmercNoBounds++; - } - // TODO: implement equatorial stereographic, see https://josm.openstreetmap.de/ticket/15970 - if (result && "stere".equals(proj) && "0".equals(parameters.get(CustomProjection.Param.lat_0.key))) { - result = false; - noEquatorStereo++; + + // exclude entries which need geodesic computation (equatorial/oblique azimuthal equidistant) + if (result && "aeqd".equals(proj)) { + final double HALF_PI = Math.PI / 2; + if (Math.abs(latitudeOfOrigin - HALF_PI) >= EPS10 && + Math.abs(latitudeOfOrigin + HALF_PI) >= EPS10) { + // See https://josm.openstreetmap.de/ticket/16129#comment:21 + result = false; + } + } + } catch (NumberFormatException | ProjectionConfigurationException e) { + e.printStackTrace(); + result = false; + } + } + + if (result && "0.0".equals(parameters.get("rf"))) { + // Proj fails with "reciprocal flattening (1/f) = 0" for + result = false; // FIXME Only for some projections? + } + + String k_0 = parameters.get("k_0"); + if (result && k_0 != null && k_0.startsWith("-")) { + // Proj fails with "k <= 0" for ESRI:102470 + result = false; } return result; } + + private static void incMap(Map map, String key) { + map.putIfAbsent(key, 0); + map.put(key, map.get(key)+1); + } } 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: + *

    + *
  • John P. Snyder (Map Projections - A Working Manual,
    + * U.S. Geological Survey Professional Paper 1395, 1987)
    • + *
  • "Coordinate Conversions and Transformations including Formulas",
    + * EPSG Guidence Note Number 7 part 2, Version 24.
    • + *
+ * + * @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; } Index: /trunk/test/unit/org/openstreetmap/josm/data/projection/ProjectionRefTest.java =================================================================== --- /trunk/test/unit/org/openstreetmap/josm/data/projection/ProjectionRefTest.java (revision 13597) +++ /trunk/test/unit/org/openstreetmap/josm/data/projection/ProjectionRefTest.java (revision 13598) @@ -79,4 +79,6 @@ static Random rand = new SecureRandom(); + static boolean debug; + /** * Setup test. @@ -92,4 +94,5 @@ */ public static void main(String[] args) throws IOException { + debug = args.length > 0 && "debug".equals(args[0]); Collection refs = readData(); refs = updateData(refs); @@ -236,14 +239,31 @@ pb.environment().put("PROJ_LIB", new File(PROJ_LIB_DIR).getAbsolutePath()); - String output; + String output = ""; try { Process process = pb.start(); OutputStream stdin = process.getOutputStream(); InputStream stdout = process.getInputStream(); + InputStream stderr = process.getErrorStream(); try (BufferedWriter writer = new BufferedWriter(new OutputStreamWriter(stdin, StandardCharsets.UTF_8))) { - writer.write(String.format("%.9f %.9f%n", ll.lon(), ll.lat())); + String s = String.format("%.9f %.9f%n", ll.lon(), ll.lat()); + if (debug) { + System.out.println("\n" + String.join(" ", args) + "\n" + s); + } + writer.write(s); } try (BufferedReader reader = new BufferedReader(new InputStreamReader(stdout, StandardCharsets.UTF_8))) { - output = reader.readLine(); + String line; + while (null != (line = reader.readLine())) { + if (debug) { + System.out.println("> " + line); + } + output = line; + } + } + try (BufferedReader reader = new BufferedReader(new InputStreamReader(stderr, StandardCharsets.UTF_8))) { + String line; + while (null != (line = reader.readLine())) { + System.err.println("! " + line); + } } } catch (IOException e) { Index: /trunk/test/unit/org/openstreetmap/josm/data/projection/ProjectionTest.java =================================================================== --- /trunk/test/unit/org/openstreetmap/josm/data/projection/ProjectionTest.java (revision 13597) +++ /trunk/test/unit/org/openstreetmap/josm/data/projection/ProjectionTest.java (revision 13598) @@ -24,4 +24,7 @@ String text; + /** + * Tests that projections are numerically stable in their definition bounds (round trip error < 1e-5) + */ @Test public void testProjections() { @@ -55,4 +58,11 @@ testProjection(Projections.getProjectionByCode("EPSG:"+Integer.toString(3942+i))); // Lambert CC9 Zones France } + + for (int i = 0; i <= 17; ++i) { + testProjection(Projections.getProjectionByCode("EPSG:"+Integer.toString(102421+i))); // WGS_1984_ARC_System Zones + } + + testProjection(Projections.getProjectionByCode("EPSG:102016")); // North Pole + testProjection(Projections.getProjectionByCode("EPSG:102019")); // South Pole if (error) { @@ -104,4 +114,7 @@ Collection projIds; + /** + * Tests that projections are numerically stable in their definition bounds (round trip error < epsilon) + */ @Test public void testProjs() { @@ -127,4 +140,7 @@ testProj("merc", 1e-5, ""); testProj("sinu", 1e-4, ""); + testProj("aeqd", 1e-5, "+lon_0=0dE +lat_0=90dN"); + testProj("aeqd", 1e-5, "+lon_0=0dE +lat_0=90dS"); + testProj("eqc", 1e-5, ""); if (error2) {