MapyCZProjection.cs


namespace GMap.NET.Projections
{
   using System;

   /// <summary>
   /// GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]
   /// PROJCS["Mapy.cz",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",134400000],PARAMETER["false_northing",-41600000],UNIT["1/32meter",0.03125]]
   /// </summary>
   public class MapyCZProjection : PureProjection
   {
      public static readonly MapyCZProjection Instance = new MapyCZProjection();

      static readonly double MinLatitude = 26;
      static readonly double MaxLatitude = 76;
      static readonly double MinLongitude = -26;
      static readonly double MaxLongitude = 38;

      #region -- Common --
      static int getLCM(int zone)
      {
         if((zone < 1) || (zone > 60))
         {
            throw new Exception("MapyCZProjection: UTM Zone number is not between 1 and 60.");
         }
         else
         {
            return ((zone * 6) - 183);
         }
      }

      static double roundoff(double xx, double yy)
      {
         var x = xx;
         var y = yy;
         x = Math.Round(x * Math.Pow(10, y)) / Math.Pow(10, y);
         return x;
      }

      static readonly double UTMSIZE = 2;
      static readonly double UNITS = 1;

      #endregion

      #region -- WGSToMapyCZ --

      public long[] WGSToPP(double la, double lo)
      {
         var utmEE = wgsToUTM(DegreesToRadians(la), DegreesToRadians(lo), 33);
         var pp = utmEEToPP(utmEE[0], utmEE[1]);
         return pp;
      }

      static long[] utmEEToPP(double east, double north)
      {
         var x = (Math.Round(east) - (-3700000.0)) * Math.Pow(2, 5);
         var y = (Math.Round(north) - (1300000.0)) * Math.Pow(2, 5);

         return new long[] { (long)x, (long)y };
      }

      double[] wgsToUTM(double la, double lo, int zone)
      {
         var latrad = la;
         var lonrad = lo;
         var latddd = RadiansToDegrees(la);
         var londdd = RadiansToDegrees(lo);

         var k = 0.9996f;
         var a = Axis;
         var f = Flattening;
         var b = a * (1.0 - f);
         var e2 = (a * a - b * b) / (a * a);
         var e = Math.Sqrt(e2);
         var ei2 = (a * a - b * b) / (b * b);
         var ei = Math.Sqrt(ei2);
         var n = (a - b) / (a + b);
         var G = a * (1.0 - n) * (1.0 - n * n) * (1.0 + (9 / 4.0) * n * n + (255.0 / 64.0) * Math.Pow(n, 4)) * (PI / 180.0);
         var w = londdd - ((double) (zone * 6 - 183));
         w = DegreesToRadians(w);
         var t = Math.Tan(latrad);
         var rho = a * (1.0 - e2) / Math.Pow(1.0 - (e2 * Math.Sin(latrad) * Math.Sin(latrad)), (3 / 2.0));
         var nu = a / Math.Sqrt(1.0 - (e2 * Math.Sin(latrad) * Math.Sin(latrad)));
         var psi = nu / rho;
         var coslat = Math.Cos(latrad);
         var sinlat = Math.Sin(latrad);
         var A0 = 1 - (e2 / 4.0) - (3 * e2 * e2 / 64.0) - (5 * Math.Pow(e2, 3) / 256.0);
         var A2 = (3 / 8.0) * (e2 + (e2 * e2 / 4.0) + (15 * Math.Pow(e2, 3) / 128.0));
         var A4 = (15 / 256.0) * (e2 * e2 + (3 * Math.Pow(e2, 3) / 4.0));
         var A6 = 35 * Math.Pow(e2, 3) / 3072.0;
         var m = a * ((A0 * latrad) - (A2 * Math.Sin(2 * latrad)) + (A4 * Math.Sin(4 * latrad)) - (A6 * Math.Sin(6 * latrad)));
         var eterm1 = (w * w / 6.0) * coslat * coslat * (psi - t * t);
         var eterm2 = (Math.Pow(w, 4) / 120.0) * Math.Pow(coslat, 4) * (4 * Math.Pow(psi, 3) * (1.0 - 6 * t * t) + psi * psi * (1.0 + 8 * t * t) - psi * 2 * t * t + Math.Pow(t, 4));
         var eterm3 = (Math.Pow(w, 6) / 5040.0) * Math.Pow(coslat, 6) * (61.0 - 479 * t * t + 179 * Math.Pow(t, 4) - Math.Pow(t, 6));
         var dE = k * nu * w * coslat * (1.0 + eterm1 + eterm2 + eterm3);
         var east = 500000.0 + (dE / UNITS);
         east = roundoff(east, UTMSIZE);
         var nterm1 = (w * w / 2.0) * nu * sinlat * coslat;
         var nterm2 = (Math.Pow(w, 4) / 24.0) * nu * sinlat * Math.Pow(coslat, 3) * (4 * psi * psi + psi - t * t);
         var nterm3 = (Math.Pow(w, 6) / 720.0) * nu * sinlat * Math.Pow(coslat, 5) * (8 * Math.Pow(psi, 4) * (11.0 - 24 * t * t) - 28 * Math.Pow(psi, 3) * (1.0 - 6 * t * t) + psi * psi * (1.0 - 32 * t * t) - psi * 2 * t * t + Math.Pow(t, 4));
         var nterm4 = (Math.Pow(w, 8) / 40320.0) * nu * sinlat * Math.Pow(coslat, 7) * (1385.0 - 3111 * t * t + 543 * Math.Pow(t, 4) - Math.Pow(t, 6));
         var dN = k * (m + nterm1 + nterm2 + nterm3 + nterm4);
         var north = (0.0 + (dN / UNITS));
         north = roundoff(north, UTMSIZE);

         return new double[] { east, north, zone };
      }

      #endregion

      #region -- MapyCZToWGS --

      public double[] PPToWGS(double x, double y)
      {
         var utmEE = ppToUTMEE(x, y);
         var ret = utmToWGS(utmEE[0], utmEE[1], 33);
         return ret;
      }

      double[] ppToUTMEE(double x, double y)
      {
         var north = y * Math.Pow(2, -5) + 1300000.0;
         var east = x * Math.Pow(2, -5) + (-3700000.0);
         east = roundoff(east, UTMSIZE);
         north = roundoff(north, UTMSIZE);

         return new double[] { east, north };
      }

      double[] utmToWGS(double eastIn, double northIn, int zone)
      {
         var k = 0.9996f;
         var a = Axis;
         var f = Flattening;
         var b = a * (1.0 - f);
         var e2 = (a * a - b * b) / (a * a);
         var e = Math.Sqrt(e2);
         var ei2 = (a * a - b * b) / (b * b);
         var ei = Math.Sqrt(ei2);
         var n = (a - b) / (a + b);
         var G = a * (1.0 - n) * (1.0 - n * n) * (1.0 + (9 / 4.0) * n * n + (255 / 64.0) * Math.Pow(n, 4)) * (PI / 180.0);
         var north = (northIn - 0) * UNITS;
         var east = (eastIn - 500000.0) * UNITS;
         var m = north / k;
         var sigma = (m * PI) / (180.0 * G);
         var footlat = sigma + ((3 * n / 2.0) - (27 * Math.Pow(n, 3) / 32.0)) * Math.Sin(2 * sigma) + ((21 * n * n / 16.0) - (55 * Math.Pow(n, 4) / 32.0)) * Math.Sin(4 * sigma) + (151 * Math.Pow(n, 3) / 96.0) * Math.Sin(6 * sigma) + (1097 * Math.Pow(n, 4) / 512.0) * Math.Sin(8 * sigma);
         var rho = a * (1.0 - e2) / Math.Pow(1.0 - (e2 * Math.Sin(footlat) * Math.Sin(footlat)), (3 / 2.0));
         var nu = a / Math.Sqrt(1.0 - (e2 * Math.Sin(footlat) * Math.Sin(footlat)));
         var psi = nu / rho;
         var t = Math.Tan(footlat);
         var x = east / (k * nu);
         var laterm1 = (t / (k * rho)) * (east * x / 2.0);
         var laterm2 = (t / (k * rho)) * (east * Math.Pow(x, 3) / 24.0) * (-4 * psi * psi + 9 * psi * (1 - t * t) + 12 * t * t);
         var laterm3 = (t / (k * rho)) * (east * Math.Pow(x, 5) / 720.0) * (8 * Math.Pow(psi, 4) * (11 - 24 * t * t) - 12 * Math.Pow(psi, 3) * (21.0 - 71 * t * t) + 15 * psi * psi * (15.0 - 98 * t * t + 15 * Math.Pow(t, 4)) + 180 * psi * (5 * t * t - 3 * Math.Pow(t, 4)) + 360 * Math.Pow(t, 4));
         var laterm4 = (t / (k * rho)) * (east * Math.Pow(x, 7) / 40320.0) * (1385.0 + 3633 * t * t + 4095 * Math.Pow(t, 4) + 1575 * Math.Pow(t, 6));
         var latrad = footlat - laterm1 + laterm2 - laterm3 + laterm4;
         var lat = RadiansToDegrees(latrad);
         var seclat = 1 / Math.Cos(footlat);
         var loterm1 = x * seclat;
         var loterm2 = (Math.Pow(x, 3) / 6.0) * seclat * (psi + 2 * t * t);
         var loterm3 = (Math.Pow(x, 5) / 120.0) * seclat * (-4 * Math.Pow(psi, 3) * (1 - 6 * t * t) + psi * psi * (9 - 68 * t * t) + 72 * psi * t * t + 24 * Math.Pow(t, 4));
         var loterm4 = (Math.Pow(x, 7) / 5040.0) * seclat * (61.0 + 662 * t * t + 1320 * Math.Pow(t, 4) + 720 * Math.Pow(t, 6));
         var w = loterm1 - loterm2 + loterm3 - loterm4;
         var longrad = DegreesToRadians(getLCM(zone)) + w;
         var lon = RadiansToDegrees(longrad);

         return new double[] { lat, lon, latrad, longrad };
      }

      #endregion

      public override RectLatLng Bounds
      {
         get
         {
            return RectLatLng.FromLTRB(MinLongitude, MaxLatitude, MaxLongitude, MinLatitude);
         }
      }

      public override GSize TileSize
      {
         get
         {
            return new GSize(256, 256);
         }
      }

      public override double Axis
      {
         get
         {
            return 6378137;
         }
      }

      public override double Flattening
      {
         get
         {
            return (1.0 / 298.257223563);
         }
      }

      public override GPoint FromLatLngToPixel(double lat, double lng, int zoom)
      {
         GPoint ret = GPoint.Empty;

         lat = Clip(lat, MinLatitude, MaxLatitude);
         lng = Clip(lng, MinLongitude, MaxLongitude);

         var size = GetTileMatrixSizePixel(zoom);
         {
            var l = WGSToPP(lat, lng);
            ret.X = (long)l[0] >> (20 - zoom);
            ret.Y = size.Height - ((long)l[1] >> (20 - zoom));
         }
         return ret;
      }

      public override PointLatLng FromPixelToLatLng(long x, long y, int zoom)
      {
         PointLatLng ret = PointLatLng.Empty;

         var size = GetTileMatrixSizePixel(zoom);

         var oX = x << (20 - zoom);
         var oY = (size.Height - y) << (20 - zoom);
         {
            var l = PPToWGS(oX, oY);
            ret.Lat = Clip(l[0], MinLatitude, MaxLatitude);
            ret.Lng = Clip(l[1], MinLongitude, MaxLongitude);
         }
         return ret;
      }

      public override GSize GetTileMatrixSizeXY(int zoom)
      {
         return new GSize((long)Math.Pow(2, zoom), (long)Math.Pow(2, zoom));
      }

      public override GSize GetTileMatrixSizePixel(int zoom)
      {
         GSize s = GetTileMatrixSizeXY(zoom);
         return new GSize(s.Width << 8, s.Height << 8);
      }

      public override GSize GetTileMatrixMinXY(int zoom)
      {
         long wh = zoom > 3 ? (3 * (long)Math.Pow(2, zoom - 4)) : 1;
         return new GSize(wh, wh);
      }

      public override GSize GetTileMatrixMaxXY(int zoom)
      {
         long wh = (long)Math.Pow(2, zoom) - (long)Math.Pow(2, zoom - 2);
         return new GSize(wh, wh);
      }
   }
}