Changeset 13598 in josm

Timestamp:

2018-04-05T19:03:04+02:00 (7 years ago)

Author:

Don-vip

Message:

see #16129 - add new projections and support for new format of ESRI file

Location:

trunk

Files:

• CONTRIBUTION (modified) (1 diff)
• scripts/BuildProjectionDefinitions.java (modified) (10 diffs)
• src/org/openstreetmap/josm/data/projection/CustomProjection.java (modified) (2 diffs)
• src/org/openstreetmap/josm/data/projection/Projections.java (modified) (5 diffs)
• src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java (modified) (1 diff)
• src/org/openstreetmap/josm/data/projection/proj/AzimuthalEquidistant.java (added)
• src/org/openstreetmap/josm/data/projection/proj/EquidistantCylindrical.java (added)
• src/org/openstreetmap/josm/data/projection/proj/ProjParameters.java (modified) (1 diff)
• test/unit/org/openstreetmap/josm/data/projection/ProjectionRefTest.java (modified) (3 diffs)
• test/unit/org/openstreetmap/josm/data/projection/ProjectionTest.java (modified) (4 diffs)

trunk/CONTRIBUTION

@@ -93 +93 @@
 Swiss CHENYX06 NTV2 grid: Source: Swiss Federal Office of Topography
 
+ESRI projection definitions: Environmental Systems Research Institute
+
 ------------------------------------ ICONS ------------------------------------
 

trunk/scripts/BuildProjectionDefinitions.java

@@ -7 +7 @@
 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 +27 @@
 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 +45 @@
     private static int noGeocent = 0;
     private static int noBaseProjection = 0;
-    private static final Map<String, Integer> baseProjectionMap = new HashMap<>();
+    private static int noEllipsoid = 0;
+    private static int noNadgrid = 0;
     private static int noDatumgrid = 0;
     private static int noJosm = 0;
@@ -47 +53 @@
     private static int noOmercNoBounds = 0;
     private static int noEquatorStereo = 0;
+
+    private static final Map<String, Integer> baseProjectionMap = new TreeMap<>();
+    private static final Map<String, Integer> ellipsoidMap = new TreeMap<>();
+    private static final Map<String, Integer> nadgridMap = new TreeMap<>();
+    private static final Map<String, Integer> datumgridMap = new TreeMap<>();
+
+    private static List<String> knownGeoidgrids;
+    private static List<String> knownNadgrids;
 
     /**
@@ -57 +71 @@
     }
 
+    static List<String> 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<String, ProjectionDefinition> 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 +87 @@
         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 +127 @@
             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 +185 @@
         }
 
-        // 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 +215 @@
         }
 
-        // 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 +237 @@
             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<String, Integer> map, String key) {
+        map.putIfAbsent(key, 0);
+        map.put(key, map.get(key)+1);
+    }
 }

trunk/src/org/openstreetmap/josm/data/projection/CustomProjection.java

@@ -296 +296 @@
             }
             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 +588 @@
             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) {

trunk/src/org/openstreetmap/josm/data/projection/Projections.java

@@ -25 +25 @@
 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 +57 @@
      */
     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 +100 @@
     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 +212 @@
     }
 
+    /**
+     * 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<? extends Proj> projClass, String origin) {
         registerBaseProjection(id, new ClassProjFactory(projClass), origin);
@@ -293 +315 @@
         List<ProjectionDefinition> 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;
     }

trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java

@@ -155 +155 @@
     }
 
+    /**
+     * 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) {

trunk/src/org/openstreetmap/josm/data/projection/proj/ProjParameters.java

@@ -2 +2 @@
 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;
 }

trunk/test/unit/org/openstreetmap/josm/data/projection/ProjectionRefTest.java

@@ -79 +79 @@
     static Random rand = new SecureRandom();
 
+    static boolean debug;
+
     /**
      * Setup test.
@@ -92 +94 @@
      */
     public static void main(String[] args) throws IOException {
+        debug = args.length > 0 && "debug".equals(args[0]);
         Collection<RefEntry> refs = readData();
         refs = updateData(refs);
@@ -236 +239 @@
         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) {

trunk/test/unit/org/openstreetmap/josm/data/projection/ProjectionTest.java

@@ -24 +24 @@
     String text;
 
+    /**
+     * Tests that projections are numerically stable in their definition bounds (round trip error < 1e-5)
+     */
     @Test
     public void testProjections() {
@@ -55 +58 @@
             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 +114 @@
     Collection<String> projIds;
 
+    /**
+     * Tests that projections are numerically stable in their definition bounds (round trip error &lt; epsilon)
+     */
     @Test
     public void testProjs() {
@@ -127 +140 @@
         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) {