TransformWgs84DutchRdNew.java

  1. package org.opentrafficsim.animation.gis;

  3. import java.awt.geom.Point2D;
  4. import java.awt.geom.Rectangle2D;

  6. /**
  7.  * Convert geographical coordinates between WGS84 and the Dutch RD (Rijksdriehoek) system. <br>
  8.  * Specific MathTransform for WGS84 (EPSG:4326) to RD_new (EPSG:28992) conversions. Code based on C code by Peter Knoppers as
  9.  * applied <a href="https://www.regiolab-delft.nl/?q=node/36">here</a>, which is based on
  10.  * <a href="https://home.solcon.nl/pvanmanen/Download/Transformatieformules.pdf">this</a> paper.
  11.  * <p>
  12.  * Copyright (c) ~2000-2024 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
  13.  * BSD-style license. See <a href="https://opentrafficsim.org/docs/license.html">OpenTrafficSim License</a>.
  14.  * </p>
  15.  * @author <a href="https://github.com/averbraeck">Alexander Verbraeck</a>
  16.  * @author <a href="https://github.com/peter-knoppers">Peter Knoppers</a>
  17.  * @author Gert-Jan Stolk
  18.  **/
  19. public final class TransformWgs84DutchRdNew
  20. {

  22.     /** Western boundary of the Dutch RD system. */
  23.     private static final double WGS84_WEST_LIMIT = 3.2;

  25.     /** Eastern boundary of the Dutch RD system. */
  26.     private static final double WGS84_EAST_LIMIT = 7.3;

  28.     /** Northern boundary of the Dutch RD system. */
  29.     private static final double WGS84_SOUTH_LIMIT = 50.6;

  31.     /** Southern boundary of the Dutch RD system. */
  32.     private static final double WGS84_NORTH_LIMIT = 53.7;

  34.     /** Western boundary of the Dutch RD system. */
  35.     private static final double RD_MINIMUM_X = 11000;

  37.     /** Eastern boundary of the Dutch RD system. */
  38.     private static final double RD_MAXIMUM_X = 280000;

  40.     /** Southern boundary of the Dutch RD system. */
  41.     private static final double RD_MINIMUM_Y = 300000;

  43.     /** Northern boundary of the Dutch RD system. */
  44.     private static final double RD_MAXIMUM_Y = 630000;

  46.     /** This class is a utility class and instances cannot be constructed. */
  47.     private TransformWgs84DutchRdNew()
  48.     {
  49.         // This class is a utility class and instances cannot be constructed.
  50.     }

  52.     /**
  53.      * Convert from WGS84 to RD coordinates.
  54.      * @param wgs84East Degrees East of Greenwich
  55.      * @param wgs84North Degrees North of the equator
  56.      * @return equivalent location in the Dutch RD system
  57.      */
  58.     private static Point2D.Double ellipseWgs84ToRd(final double wgs84East, final double wgs84North)
  59.     {
  60.         if (wgs84North > WGS84_NORTH_LIMIT || wgs84North < WGS84_SOUTH_LIMIT || wgs84East < WGS84_WEST_LIMIT
  61.                 || (wgs84East > WGS84_EAST_LIMIT))
  62.         {
  63.             throw new IllegalArgumentException("ellipswgs842rd input out of range (" + wgs84East + ", " + wgs84North + ")");
  64.         }
  65.         //@formatter:off
  66.         @SuppressWarnings("checkstyle:nowhitespacebefore")
  67.         /** Coefficients for ellipsewgs842rd. */
  68.         final double[][] r = 
  69.         { /* p down, q right */
  70.             {  155000.00, 190094.945,   -0.008, -32.391, 0.0   , },
  71.             {     -0.705, -11832.228,    0.0  ,   0.608, 0.0   , },
  72.             {      0.0  ,   -114.221,    0.0  ,   0.148, 0.0   , },
  73.             {      0.0  ,     -2.340,    0.0  ,   0.0  , 0.0   , },
  74.             {      0.0  ,      0.0  ,    0.0  ,   0.0  , 0.0   , }
  75.         };
  76.         @SuppressWarnings("checkstyle:nowhitespacebefore")
  77.         final double[][] s = 
  78.         { /* p down, q right */
  79.             { 463000.00 ,      0.433, 3638.893,   0.0  ,  0.092, },
  80.             { 309056.544,     -0.032, -157.984,   0.0  , -0.054, },
  81.             {     73.077,      0.0  ,   -6.439,   0.0  ,  0.0  , },
  82.             {     59.788,      0.0  ,    0.0  ,   0.0  ,  0.0  , },
  83.             {      0.0  ,      0.0  ,    0.0  ,   0.0  ,  0.0  , }
  84.         };
  85.         //@formatter:on
  86.         double resultX = 0;
  87.         double resultY = 0;
  88.         double powNorth = 1;
  89.         double dNorth = 0.36 * (wgs84North - 52.15517440);
  90.         double dEast = 0.36 * (wgs84East - 5.38720621);

  92.         for (int p = 0; p < 5; p++)
  93.         {
  94.             double powEast = 1;
  95.             for (int q = 0; q < 5; q++)
  96.             {
  97.                 resultX += r[p][q] * powEast * powNorth;
  98.                 resultY += s[p][q] * powEast * powNorth;
  99.                 powEast *= dEast;
  100.             }
  101.             powNorth *= dNorth;
  102.         }
  103.         return new Point2D.Double(resultX, resultY);
  104.     }

  106.     /**
  107.      * Convert coordinates from WGS84 to the Dutch RD system.
  108.      * @param rdX X coordinate in the Dutch RD system
  109.      * @param rdY Y coordinate in the Dutch RD system
  110.      * @return equivalent location in the WGS84 system
  111.      * @throws IllegalArgumentException when the given coordinates are not within the area of the Dutch RD system
  112.      */
  113.     private static Point2D rdToEllipseWgs84(final double rdX, final double rdY) throws IllegalArgumentException
  114.     {
  115.         if (rdX < RD_MINIMUM_X || rdX > RD_MAXIMUM_X || rdY < RD_MINIMUM_Y || rdY > RD_MAXIMUM_Y)
  116.         {
  117.             throw new IllegalArgumentException(
  118.                     "Location (" + rdX + "," + rdY + ") is not within the range " + "of the Dutch RD system");
  119.         }
  120.         final double dX = (rdX - 155000) / 100000;
  121.         final double dY = (rdY - 463000) / 100000;
  122.         /* Coefficients are for zone 31 (0E .. 6E); roughly west of Apeldoorn */
  123.         //@formatter:off
  124.         @SuppressWarnings("checkstyle:nowhitespacebefore")
  125.         final double[][] k = 
  126.         { /* p down, q right */
  127.             { 3600 * 52.15517440, 3235.65389, -0.24750, -0.06550, 0.0    , },
  128.             {        -0.00738   ,   -0.00012,  0.0    ,  0.0    , 0.0    , },
  129.             {       -32.58297   ,   -0.84978, -0.01709, -0.00039, 0.0    , },
  130.             {         0.0       ,    0.0    ,  0.0    ,  0.0    , 0.0    , },
  131.             {         0.00530   ,    0.00033,  0.0    ,  0.0    , 0.0    , },
  132.             {         0.0       ,    0.0    ,  0.0    ,  0.0    , 0.0    , }
  133.         };
  134.         @SuppressWarnings("checkstyle:nowhitespacebefore")
  135.         final double[][] l = { /* p down, q right */
  136.             {  3600 * 5.38720621,    0.01199,  0.00022,  0.0    , 0.0    , },
  137.             {      5260.52916   ,  105.94684,  2.45656,  0.05594, 0.00128, },
  138.             {        -0.00022   ,    0.0    ,  0.0    ,  0.0    , 0.0    , },
  139.             {        -0.81885   ,   -0.05607, -0.00256,  0.0    , 0.0    , },
  140.             {         0.0       ,    0.0    ,  0.0    ,  0.0    , 0.0    , },
  141.             {         0.00026   ,    0.0    ,  0.0    ,  0.0    , 0.0    , }
  142.         };
  143.         //@formatter:on
  144.         double resultNorth = 0;
  145.         double resultEast = 0;
  146.         double powX = 1;
  147.         for (int p = 0; p < 6; p++)
  148.         {
  149.             double powY = 1;
  150.             for (int q = 0; q < 5; q++)
  151.             {
  152.                 resultNorth += k[p][q] * powX * powY / 3600;
  153.                 resultEast += l[p][q] * powX * powY / 3600;
  154.                 powY *= dY;
  155.             }
  156.             powX *= dX;
  157.         }
  158.         return new Point2D.Double(resultEast, resultNorth);
  159.     }

  161.     /**
  162.      * Convert a coordinate pair in the local system to WGS84 coordinates.
  163.      * @param local coordinates in the local system.
  164.      * @return the equivalent location in degrees in the WGS84 coordinate system
  165.      * @throws IllegalArgumentException when <cite>local</cite> is not valid in the local system
  166.      */
  167.     public static Point2D toWgs84(final Point2D local) throws IllegalArgumentException
  168.     {
  169.         return toWgs84(local.getX(), local.getY());
  170.     }

  172.     /**
  173.      * Convert a coordinate pair in the local system to WGS84 coordinates.
  174.      * @param localX X-coordinate in the local system
  175.      * @param localY Y-coordinate in the local system
  176.      * @return the equivalent location in degrees in the WGS84 coordinate system
  177.      * @throws IllegalArgumentException when <cite>localX</cite>, <cite>localY</cite> is not valid in the local system
  178.      */
  179.     public static Point2D toWgs84(final double localX, final double localY) throws IllegalArgumentException
  180.     {
  181.         return rdToEllipseWgs84(localX, localY);
  182.     }

  184.     /**
  185.      * Convert a coordinate pair in WGS84 coordinates to local coordinates.
  186.      * @param wgs84 coordinates in degrees in the WGS84 coordinate system
  187.      * @return the equivalent location in the local coordinate system
  188.      * @throws IllegalArgumentException when <cite>wgs84</cite> is not valid in the local system
  189.      */
  190.     public static Point2D fromWgs84(final Point2D wgs84) throws IllegalArgumentException
  191.     {
  192.         return fromWgs84(wgs84.getX(), wgs84.getY());
  193.     }

  195.     /**
  196.      * Convert a coordinate pair in WGS84 coordinates to local coordinates.
  197.      * @param wgs84East East coordinate in degrees in the WGS84 system (negative value indicates West)
  198.      * @param wgs84North North coordinate in degrees in the WGS84 system (negative value indicates South)
  199.      * @return the equivalent location in the local coordinate system
  200.      * @throws IllegalArgumentException when <cite>wgs84</cite> is not valid in the local system
  201.      */
  202.     public static Point2D fromWgs84(final double wgs84East, final double wgs84North) throws IllegalArgumentException
  203.     {
  204.         return ellipseWgs84ToRd(wgs84East, wgs84North);
  205.     }

  207.     /**
  208.      * Report the bounding box for conversion to the local coordinate system.<br>
  209.      * Conversions from WGS84 to the local coordinate system should fail for locations outside this bounding box. If the valid
  210.      * range is not adequately described by a rectangular bounding box, conversions for some areas within this bounding box may
  211.      * also fail (with an IllegalArgumentException). There is no guarantee that the result of a conversion lies within the
  212.      * bounding box for the reverse conversion.
  213.      * @return bounding box in WGS84 degrees
  214.      */
  215.     public static Rectangle2D fromWgs84Bounds()
  216.     {
  217.         return new Rectangle2D.Double(WGS84_WEST_LIMIT, WGS84_SOUTH_LIMIT, WGS84_EAST_LIMIT - WGS84_WEST_LIMIT,
  218.                 WGS84_NORTH_LIMIT - WGS84_SOUTH_LIMIT);
  219.     }

  221.     /**
  222.      * Report the bounding box for conversions from the local coordinate system. <br>
  223.      * Conversions from the local coordinate system to WGS84 should fail for locations outside this bounding box. If the valid
  224.      * range is not adequately described by a rectangular bounding box, conversions for some areas within this bounding box may
  225.      * also fail (with an IllegalArgumentException). There is no guarantee that the result of a conversion lies within the
  226.      * bounding box for the reverse conversion.
  227.      * @return bounding box of the local coordinate system
  228.      */
  229.     public static Rectangle2D toWgs84Bounds()
  230.     {
  231.         return new Rectangle2D.Double(RD_MINIMUM_X, RD_MINIMUM_Y, RD_MAXIMUM_X - RD_MINIMUM_X, RD_MAXIMUM_Y - RD_MINIMUM_Y);
  232.     }

  234. }
Created with JaCoCo 0.8.12.202403310830