1 /* 2 * This file is part of the "GKMap". 3 * GKMap project borrowed from GMap.NET (by radioman). 4 * 5 * Copyright (C) 2009-2018 by radioman (email@radioman.lt). 6 * This program is licensed under the FLAT EARTH License. 7 */ 8 9 using System; 10 using System.Collections.Generic; 11 using System.Diagnostics; 12 13 namespace GKMap 14 { 15 /// <summary> 16 /// defines projection 17 /// </summary> 18 public abstract class PureProjection 19 { 20 private readonly List<Dictionary<PointLatLng, GPoint>> FromLatLngToPixelCache = new List<Dictionary<PointLatLng, GPoint>>(33); 21 private readonly List<Dictionary<GPoint, PointLatLng>> FromPixelToLatLngCache = new List<Dictionary<GPoint, PointLatLng>>(33); 22 23 /// <summary> 24 /// size of tile 25 /// </summary> 26 public abstract GSize TileSize 27 { 28 get; 29 } 30 31 /// <summary> 32 /// Semi-major axis of ellipsoid, in meters 33 /// </summary> 34 public abstract double Axis 35 { 36 get; 37 } 38 39 /// <summary> 40 /// Flattening of ellipsoid 41 /// </summary> 42 public abstract double Flattening 43 { 44 get; 45 } 46 47 PureProjection()48 protected PureProjection() 49 { 50 for (int i = 0; i < FromLatLngToPixelCache.Capacity; i++) { 51 FromLatLngToPixelCache.Add(new Dictionary<PointLatLng, GPoint>()); 52 FromPixelToLatLngCache.Add(new Dictionary<GPoint, PointLatLng>()); 53 } 54 } 55 56 /// <summary> 57 /// get pixel coordinates from lat/lng 58 /// </summary> 59 /// <param name="lat"></param> 60 /// <param name="lng"></param> 61 /// <param name="zoom"></param> 62 /// <returns></returns> FromLatLngToPixel(double lat, double lng, int zoom)63 public abstract GPoint FromLatLngToPixel(double lat, double lng, int zoom); 64 65 /// <summary> 66 /// gets lat/lng coordinates from pixel coordinates 67 /// </summary> 68 /// <param name="x"></param> 69 /// <param name="y"></param> 70 /// <param name="zoom"></param> 71 /// <returns></returns> FromPixelToLatLng(long x, long y, int zoom)72 public abstract PointLatLng FromPixelToLatLng(long x, long y, int zoom); 73 74 /// <summary> 75 /// get pixel coordinates from lat/lng 76 /// </summary> 77 /// <param name="p"></param> 78 /// <param name="zoom"></param> 79 /// <returns></returns> FromLatLngToPixel(PointLatLng p, int zoom, bool useCache = false)80 public GPoint FromLatLngToPixel(PointLatLng p, int zoom, bool useCache = false) 81 { 82 if (useCache) { 83 GPoint ret; 84 if (!FromLatLngToPixelCache[zoom].TryGetValue(p, out ret)) { 85 ret = FromLatLngToPixel(p.Lat, p.Lng, zoom); 86 FromLatLngToPixelCache[zoom].Add(p, ret); 87 88 // for reverse cache 89 if (!FromPixelToLatLngCache[zoom].ContainsKey(ret)) { 90 FromPixelToLatLngCache[zoom].Add(ret, p); 91 } 92 93 Debug.WriteLine("FromLatLngToPixelCache[" + zoom + "] added " + p + " with " + ret); 94 } 95 return ret; 96 } else { 97 return FromLatLngToPixel(p.Lat, p.Lng, zoom); 98 } 99 } 100 101 /// <summary> 102 /// gets lat/lng coordinates from pixel coordinates 103 /// </summary> 104 /// <param name="p"></param> 105 /// <param name="zoom"></param> 106 /// <returns></returns> FromPixelToLatLng(GPoint p, int zoom, bool useCache = false)107 public PointLatLng FromPixelToLatLng(GPoint p, int zoom, bool useCache = false) 108 { 109 if (useCache) { 110 PointLatLng ret; 111 if (!FromPixelToLatLngCache[zoom].TryGetValue(p, out ret)) { 112 ret = FromPixelToLatLng(p.X, p.Y, zoom); 113 FromPixelToLatLngCache[zoom].Add(p, ret); 114 115 // for reverse cache 116 if (!FromLatLngToPixelCache[zoom].ContainsKey(ret)) { 117 FromLatLngToPixelCache[zoom].Add(ret, p); 118 } 119 120 Debug.WriteLine("FromPixelToLatLngCache[" + zoom + "] added " + p + " with " + ret); 121 } 122 return ret; 123 } else { 124 return FromPixelToLatLng(p.X, p.Y, zoom); 125 } 126 } 127 128 /// <summary> 129 /// gets tile coordinate from pixel coordinates 130 /// </summary> 131 /// <param name="p"></param> 132 /// <returns></returns> FromPixelToTileXY(GPoint p)133 public virtual GPoint FromPixelToTileXY(GPoint p) 134 { 135 return new GPoint(p.X / TileSize.Width, p.Y / TileSize.Height); 136 } 137 138 /// <summary> 139 /// gets pixel coordinate from tile coordinate 140 /// </summary> 141 /// <param name="p"></param> 142 /// <returns></returns> FromTileXYToPixel(GPoint p)143 public virtual GPoint FromTileXYToPixel(GPoint p) 144 { 145 return new GPoint((p.X * TileSize.Width), (p.Y * TileSize.Height)); 146 } 147 148 /// <summary> 149 /// min. tile in tiles at custom zoom level 150 /// </summary> 151 /// <param name="zoom"></param> 152 /// <returns></returns> GetTileMatrixMinXY(int zoom)153 public abstract GSize GetTileMatrixMinXY(int zoom); 154 155 /// <summary> 156 /// max. tile in tiles at custom zoom level 157 /// </summary> 158 /// <param name="zoom"></param> 159 /// <returns></returns> GetTileMatrixMaxXY(int zoom)160 public abstract GSize GetTileMatrixMaxXY(int zoom); 161 162 /// <summary> 163 /// gets matrix size in tiles 164 /// </summary> 165 /// <param name="zoom"></param> 166 /// <returns></returns> GetTileMatrixSizeXY(int zoom)167 public virtual GSize GetTileMatrixSizeXY(int zoom) 168 { 169 GSize sMin = GetTileMatrixMinXY(zoom); 170 GSize sMax = GetTileMatrixMaxXY(zoom); 171 172 return new GSize(sMax.Width - sMin.Width + 1, sMax.Height - sMin.Height + 1); 173 } 174 175 /// <summary> 176 /// tile matrix size in pixels at custom zoom level 177 /// </summary> 178 /// <param name="zoom"></param> 179 /// <returns></returns> GetTileMatrixItemCount(int zoom)180 public long GetTileMatrixItemCount(int zoom) 181 { 182 GSize s = GetTileMatrixSizeXY(zoom); 183 return (s.Width * s.Height); 184 } 185 186 /// <summary> 187 /// gets matrix size in pixels 188 /// </summary> 189 /// <param name="zoom"></param> 190 /// <returns></returns> GetTileMatrixSizePixel(int zoom)191 public virtual GSize GetTileMatrixSizePixel(int zoom) 192 { 193 GSize s = GetTileMatrixSizeXY(zoom); 194 return new GSize(s.Width * TileSize.Width, s.Height * TileSize.Height); 195 } 196 197 /// <summary> 198 /// gets all tiles in rect at specific zoom 199 /// </summary> GetAreaTileList(RectLatLng rect, int zoom, int padding)200 public List<GPoint> GetAreaTileList(RectLatLng rect, int zoom, int padding) 201 { 202 List<GPoint> ret = new List<GPoint>(); 203 204 GPoint topLeft = FromPixelToTileXY(FromLatLngToPixel(rect.LocationTopLeft, zoom)); 205 GPoint rightBottom = FromPixelToTileXY(FromLatLngToPixel(rect.LocationRightBottom, zoom)); 206 207 for (long x = (topLeft.X - padding); x <= (rightBottom.X + padding); x++) { 208 for (long y = (topLeft.Y - padding); y <= (rightBottom.Y + padding); y++) { 209 GPoint p = new GPoint(x, y); 210 if (!ret.Contains(p) && p.X >= 0 && p.Y >= 0) { 211 ret.Add(p); 212 } 213 } 214 } 215 216 return ret; 217 } 218 219 /// <summary> 220 /// The ground resolution indicates the distance (in meters) on the ground that’s represented by a single pixel in the map. 221 /// For example, at a ground resolution of 10 meters/pixel, each pixel represents a ground distance of 10 meters. 222 /// </summary> 223 /// <param name="zoom"></param> 224 /// <param name="latitude"></param> 225 /// <returns></returns> GetGroundResolution(int zoom, double latitude)226 public virtual double GetGroundResolution(int zoom, double latitude) 227 { 228 return (Math.Cos(latitude * (Math.PI / 180)) * 2 * Math.PI * Axis) / GetTileMatrixSizePixel(zoom).Width; 229 } 230 231 /// <summary> 232 /// gets boundaries 233 /// </summary> 234 public virtual RectLatLng Bounds 235 { 236 get { 237 return RectLatLng.FromLTRB(-180, 90, 180, -90); 238 } 239 } 240 241 #region -- math functions -- 242 243 /// <summary> 244 /// PI 245 /// </summary> 246 protected static readonly double PI = Math.PI; 247 248 /// <summary> 249 /// Half of PI 250 /// </summary> 251 protected static readonly double HALF_PI = (PI * 0.5); 252 253 /// <summary> 254 /// PI * 2 255 /// </summary> 256 protected static readonly double TWO_PI = (PI * 2.0); 257 258 /// <summary> 259 /// EPSLoN 260 /// </summary> 261 protected static readonly double EPSLoN = 1.0e-10; 262 263 /// <summary> 264 /// MAX_VAL 265 /// </summary> 266 protected const double MAX_VAL = 4; 267 268 /// <summary> 269 /// MAXLONG 270 /// </summary> 271 protected static readonly double MAXLONG = 2147483647; 272 273 /// <summary> 274 /// DBLLONG 275 /// </summary> 276 protected static readonly double DBLLONG = 4.61168601e18; 277 278 static readonly double R2D = 180 / Math.PI; 279 static readonly double D2R = Math.PI / 180; 280 DegreesToRadians(double deg)281 public static double DegreesToRadians(double deg) 282 { 283 return (D2R * deg); 284 } 285 RadiansToDegrees(double rad)286 public static double RadiansToDegrees(double rad) 287 { 288 return (R2D * rad); 289 } 290 291 ///<summary> 292 /// return the sign of an argument 293 ///</summary> Sign(double x)294 protected static double Sign(double x) 295 { 296 if (x < 0.0) 297 return (-1); 298 else 299 return (1); 300 } 301 AdjustLongitude(double x)302 protected static double AdjustLongitude(double x) 303 { 304 long count = 0; 305 while (true) { 306 if (Math.Abs(x) <= PI) 307 break; 308 else 309 if (((long)Math.Abs(x / Math.PI)) < 2) 310 x = x - (Sign(x) * TWO_PI); 311 else 312 if (((long)Math.Abs(x / TWO_PI)) < MAXLONG) { 313 x = x - (((long)(x / TWO_PI)) * TWO_PI); 314 } else 315 if (((long)Math.Abs(x / (MAXLONG * TWO_PI))) < MAXLONG) { 316 x = x - (((long)(x / (MAXLONG * TWO_PI))) * (TWO_PI * MAXLONG)); 317 } else 318 if (((long)Math.Abs(x / (DBLLONG * TWO_PI))) < MAXLONG) { 319 x = x - (((long)(x / (DBLLONG * TWO_PI))) * (TWO_PI * DBLLONG)); 320 } else 321 x = x - (Sign(x) * TWO_PI); 322 count++; 323 if (count > MAX_VAL) 324 break; 325 } 326 return (x); 327 } 328 329 /// <summary> 330 /// calculates the sine and cosine 331 /// </summary> SinCos(double val, out double sin, out double cos)332 protected static void SinCos(double val, out double sin, out double cos) 333 { 334 sin = Math.Sin(val); 335 cos = Math.Cos(val); 336 } 337 338 /// <summary> 339 /// computes the constants e0, e1, e2, and e3 which are used 340 /// in a series for calculating the distance along a meridian. 341 /// </summary> 342 /// <param name="x">represents the eccentricity squared</param> 343 /// <returns></returns> e0fn(double x)344 protected static double e0fn(double x) 345 { 346 return (1.0 - 0.25 * x * (1.0 + x / 16.0 * (3.0 + 1.25 * x))); 347 } 348 e1fn(double x)349 protected static double e1fn(double x) 350 { 351 return (0.375 * x * (1.0 + 0.25 * x * (1.0 + 0.46875 * x))); 352 } 353 e2fn(double x)354 protected static double e2fn(double x) 355 { 356 return (0.05859375 * x * x * (1.0 + 0.75 * x)); 357 } 358 e3fn(double x)359 protected static double e3fn(double x) 360 { 361 return (x * x * x * (35.0 / 3072.0)); 362 } 363 364 /// <summary> 365 /// computes the value of M which is the distance along a meridian 366 /// from the Equator to latitude phi. 367 /// </summary> mlfn(double e0, double e1, double e2, double e3, double phi)368 protected static double mlfn(double e0, double e1, double e2, double e3, double phi) 369 { 370 return (e0 * phi - e1 * Math.Sin(2.0 * phi) + e2 * Math.Sin(4.0 * phi) - e3 * Math.Sin(6.0 * phi)); 371 } 372 373 /// <summary> 374 /// calculates UTM zone number 375 /// </summary> 376 /// <param name="lon">Longitude in degrees</param> 377 /// <returns></returns> GetUTMzone(double lon)378 protected static long GetUTMzone(double lon) 379 { 380 return ((long)(((lon + 180.0) / 6.0) + 1.0)); 381 } 382 383 /// <summary> 384 /// Clips a number to the specified minimum and maximum values. 385 /// </summary> 386 /// <param name="n">The number to clip.</param> 387 /// <param name="minValue">Minimum allowable value.</param> 388 /// <param name="maxValue">Maximum allowable value.</param> 389 /// <returns>The clipped value.</returns> Clip(double n, double minValue, double maxValue)390 protected static double Clip(double n, double minValue, double maxValue) 391 { 392 return Math.Min(Math.Max(n, minValue), maxValue); 393 } 394 395 /// <summary> 396 /// distance (in km) between two points specified by latitude/longitude 397 /// The Haversine formula, http://www.movable-type.co.uk/scripts/latlong.html 398 /// </summary> 399 /// <param name="p1"></param> 400 /// <param name="p2"></param> 401 /// <returns></returns> GetDistance(PointLatLng p1, PointLatLng p2)402 public double GetDistance(PointLatLng p1, PointLatLng p2) 403 { 404 double dLat1InRad = p1.Lat * (Math.PI / 180); 405 double dLong1InRad = p1.Lng * (Math.PI / 180); 406 double dLat2InRad = p2.Lat * (Math.PI / 180); 407 double dLong2InRad = p2.Lng * (Math.PI / 180); 408 double dLongitude = dLong2InRad - dLong1InRad; 409 double dLatitude = dLat2InRad - dLat1InRad; 410 double a = Math.Pow(Math.Sin(dLatitude / 2), 2) + Math.Cos(dLat1InRad) * Math.Cos(dLat2InRad) * Math.Pow(Math.Sin(dLongitude / 2), 2); 411 double c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a)); 412 double dDistance = (Axis / 1000.0) * c; 413 return dDistance; 414 } 415 GetDistanceInPixels(GPoint point1, GPoint point2)416 public double GetDistanceInPixels(GPoint point1, GPoint point2) 417 { 418 double a = point2.X - point1.X; 419 double b = point2.Y - point1.Y; 420 421 return Math.Sqrt(a * a + b * b); 422 } 423 424 /// <summary> 425 /// Accepts two coordinates in degrees. 426 /// </summary> 427 /// <returns>A double value in degrees. From 0 to 360.</returns> GetBearing(PointLatLng p1, PointLatLng p2)428 public double GetBearing(PointLatLng p1, PointLatLng p2) 429 { 430 var latitude1 = DegreesToRadians(p1.Lat); 431 var latitude2 = DegreesToRadians(p2.Lat); 432 var longitudeDifference = DegreesToRadians(p2.Lng - p1.Lng); 433 434 var y = Math.Sin(longitudeDifference) * Math.Cos(latitude2); 435 var x = Math.Cos(latitude1) * Math.Sin(latitude2) - Math.Sin(latitude1) * Math.Cos(latitude2) * Math.Cos(longitudeDifference); 436 437 return (RadiansToDegrees(Math.Atan2(y, x)) + 360) % 360; 438 } 439 440 /// <summary> 441 /// Conversion from cartesian earth-sentered coordinates to geodetic coordinates in the given datum 442 /// </summary> 443 /// <param name="Lat"></param> 444 /// <param name="Lon"></param> 445 /// <param name="Height">Height above ellipsoid [m]</param> 446 /// <param name="X"></param> 447 /// <param name="Y"></param> 448 /// <param name="Z"></param> FromGeodeticToCartesian(double Lat, double Lng, double Height, out double X, out double Y, out double Z)449 public void FromGeodeticToCartesian(double Lat, double Lng, double Height, out double X, out double Y, out double Z) 450 { 451 Lat = (Math.PI / 180) * Lat; 452 Lng = (Math.PI / 180) * Lng; 453 454 double B = Axis * (1.0 - Flattening); 455 double ee = 1.0 - (B / Axis) * (B / Axis); 456 double N = (Axis / Math.Sqrt(1.0 - ee * Math.Sin(Lat) * Math.Sin(Lat))); 457 458 X = (N + Height) * Math.Cos(Lat) * Math.Cos(Lng); 459 Y = (N + Height) * Math.Cos(Lat) * Math.Sin(Lng); 460 Z = (N * (B / Axis) * (B / Axis) + Height) * Math.Sin(Lat); 461 } 462 463 /// <summary> 464 /// Conversion from cartesian earth-sentered coordinates to geodetic coordinates in the given datum 465 /// </summary> 466 /// <param name="X"></param> 467 /// <param name="Y"></param> 468 /// <param name="Z"></param> 469 /// <param name="Lat"></param> 470 /// <param name="Lon"></param> FromCartesianTGeodetic(double X, double Y, double Z, out double Lat, out double Lng)471 public void FromCartesianTGeodetic(double X, double Y, double Z, out double Lat, out double Lng) 472 { 473 double E = Flattening * (2.0 - Flattening); 474 Lng = Math.Atan2(Y, X); 475 476 double P = Math.Sqrt(X * X + Y * Y); 477 double Theta = Math.Atan2(Z, (P * (1.0 - Flattening))); 478 double st = Math.Sin(Theta); 479 double ct = Math.Cos(Theta); 480 Lat = Math.Atan2(Z + E / (1.0 - Flattening) * Axis * st * st * st, P - E * Axis * ct * ct * ct); 481 482 Lat /= (Math.PI / 180); 483 Lng /= (Math.PI / 180); 484 } 485 486 #endregion 487 } 488 } 489