地理座標變換

Question

我想一個類，可以在一個系統轉換爲另一種。

我發現Python中的源代碼，並試圖將它移植到C＃。

這是Python源。 From here

import math 

class GlobalMercator(object): 

    def __init__(self, tileSize=256): 
     "Initialize the TMS Global Mercator pyramid" 
     self.tileSize = tileSize 
     self.initialResolution = 2 * math.pi * 6378137/self.tileSize 
     # 156543.03392804062 for tileSize 256 pixels 
     self.originShift = 2 * math.pi * 6378137/2.0 
     # 20037508.342789244 

    def LatLonToMeters(self, lat, lon): 
     "Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913" 

     mx = lon * self.originShift/180.0 
     my = math.log(math.tan((90 + lat) * math.pi/360.0))/(math.pi/180.0) 

     my = my * self.originShift/180.0 
     return mx, my 

    def MetersToLatLon(self, mx, my): 
     "Converts XY point from Spherical Mercator EPSG:900913 to lat/lon in WGS84 Datum" 

     lon = (mx/self.originShift) * 180.0 
     lat = (my/self.originShift) * 180.0 

     lat = 180/math.pi * (2 * math.atan(math.exp(lat * math.pi/180.0)) - math.pi/2.0) 
     return lat, lon 

    def PixelsToMeters(self, px, py, zoom): 
     "Converts pixel coordinates in given zoom level of pyramid to EPSG:900913" 

     res = self.Resolution(zoom) 
     mx = px * res - self.originShift 
     my = py * res - self.originShift 
     return mx, my 

    def MetersToPixels(self, mx, my, zoom): 
     "Converts EPSG:900913 to pyramid pixel coordinates in given zoom level" 

     res = self.Resolution(zoom) 
     px = (mx + self.originShift)/res 
     py = (my + self.originShift)/res 
     return px, py 

    def PixelsToTile(self, px, py): 
     "Returns a tile covering region in given pixel coordinates" 

     tx = int(math.ceil(px/float(self.tileSize)) - 1) 
     ty = int(math.ceil(py/float(self.tileSize)) - 1) 
     return tx, ty 

    def PixelsToRaster(self, px, py, zoom): 
     "Move the origin of pixel coordinates to top-left corner" 

     mapSize = self.tileSize << zoom 
     return px, mapSize - py 

    def MetersToTile(self, mx, my, zoom): 
     "Returns tile for given mercator coordinates" 

     px, py = self.MetersToPixels(mx, my, zoom) 
     return self.PixelsToTile(px, py) 

    def TileBounds(self, tx, ty, zoom): 
     "Returns bounds of the given tile in EPSG:900913 coordinates" 

     minx, miny = self.PixelsToMeters(tx*self.tileSize, ty*self.tileSize, zoom) 
     maxx, maxy = self.PixelsToMeters((tx+1)*self.tileSize, (ty+1)*self.tileSize, zoom) 
     return (minx, miny, maxx, maxy) 

    def TileLatLonBounds(self, tx, ty, zoom): 
     "Returns bounds of the given tile in latutude/longitude using WGS84 datum" 

     bounds = self.TileBounds(tx, ty, zoom) 
     minLat, minLon = self.MetersToLatLon(bounds[0], bounds[1]) 
     maxLat, maxLon = self.MetersToLatLon(bounds[2], bounds[3]) 

     return (minLat, minLon, maxLat, maxLon) 

    def Resolution(self, zoom): 
     "Resolution (meters/pixel) for given zoom level (measured at Equator)" 

     # return (2 * math.pi * 6378137)/(self.tileSize * 2**zoom) 
     return self.initialResolution/(2**zoom) 

    def ZoomForPixelSize(self, pixelSize): 
     "Maximal scaledown zoom of the pyramid closest to the pixelSize." 

     for i in range(30): 
      if pixelSize > self.Resolution(i): 
       return i-1 if i!=0 else 0 # We don't want to scale up 

    def GoogleTile(self, tx, ty, zoom): 
     "Converts TMS tile coordinates to Google Tile coordinates" 

     # coordinate origin is moved from bottom-left to top-left corner of the extent 
     return tx, (2**zoom - 1) - ty 

    def QuadTree(self, tx, ty, zoom): 
     "Converts TMS tile coordinates to Microsoft QuadTree" 

     quadKey = "" 
     ty = (2**zoom - 1) - ty 
     for i in range(zoom, 0, -1): 
      digit = 0 
      mask = 1 << (i-1) 
      if (tx & mask) != 0: 
       digit += 1 
      if (ty & mask) != 0: 
       digit += 2 
      quadKey += str(digit) 

     return quadKey 

這是我的C＃端口。

public class GlobalMercator { 
    private Int32 TileSize; 
    private Double InitialResolution; 
    private Double OriginShift; 
    private const Int32 EarthRadius = 6378137; 
    public GlobalMercator() { 
     TileSize = 256; 
     InitialResolution = 2 * Math.PI * EarthRadius/TileSize; 
     OriginShift = 2 * Math.PI * EarthRadius/2; 
    } 
    public DPoint LatLonToMeters(Double lat, Double lon) { 
     var p = new DPoint(); 
     p.X = lon * OriginShift/180; 
     p.Y = Math.Log(Math.Tan((90 + lat) * Math.PI/360))/(Math.PI/180); 
     p.Y = p.Y * OriginShift/180; 
     return p; 
    } 
    public GeoPoint MetersToLatLon(DPoint m) { 
     var ll = new GeoPoint(); 
     ll.Longitude = (m.X/OriginShift) * 180; 
     ll.Latitude = (m.Y/OriginShift) * 180; 
     ll.Latitude = 180/Math.PI * (2 * Math.Atan(Math.Exp(ll.Latitude * Math.PI/180)) - Math.PI/2); 
     return ll; 
    } 
    public DPoint PixelsToMeters(DPoint p, Int32 zoom) { 
     var res = Resolution(zoom); 
     var met = new DPoint(); 
     met.X = p.X * res - OriginShift; 
     met.Y = p.Y * res - OriginShift; 
     return met; 
    } 
    public DPoint MetersToPixels(DPoint m, Int32 zoom) { 
     var res = Resolution(zoom); 
     var pix = new DPoint(); 
     pix.X = (m.X + OriginShift)/res; 
     pix.Y = (m.Y + OriginShift)/res; 
     return pix; 
    } 
    public Point PixelsToTile(DPoint p) { 
     var t = new Point(); 
     t.X = (Int32)Math.Ceiling(p.X/(Double)TileSize) - 1; 
     t.Y = (Int32)Math.Ceiling(p.Y/(Double)TileSize) - 1; 
     return t; 
    } 
    public Point PixelsToRaster(Point p, Int32 zoom) { 
     var mapSize = TileSize << zoom; 
     return new Point(p.X, mapSize - p.Y); 
    } 
    public Point MetersToTile(Point m, Int32 zoom) { 
     var p = MetersToPixels(m, zoom); 
     return PixelsToTile(p); 
    } 

    public Pair<DPoint> TileBounds(Point t, Int32 zoom) { 
     var min = PixelsToMeters(new DPoint(t.X * TileSize, t.Y * TileSize), zoom); 
     var max = PixelsToMeters(new DPoint((t.X + 1) * TileSize, (t.Y + 1) * TileSize), zoom); 
     return new Pair<DPoint>(min, max); 
    } 
    public Pair<GeoPoint> TileLatLonBounds(Point t, Int32 zoom) { 
     var bound = TileBounds(t, zoom); 
     var min = MetersToLatLon(bound.Min); 
     var max = MetersToLatLon(bound.Max); 
     return new Pair<GeoPoint>(min, max); 
    } 
    public Double Resolution(Int32 zoom) { 
     return InitialResolution/(2^zoom); 
    } 
    public Double ZoomForPixelSize(Double pixelSize) { 
     for (var i = 0; i < 30; i++) 
      if (pixelSize > Resolution(i)) 
       return i != 0 ? i - 1 : 0; 
     throw new InvalidOperationException(); 
    } 
    public Point ToGoogleTile(Point t, Int32 zoom) { 
     return new Point(t.X, ((Int32)Math.Pow(2, zoom) - 1) - t.Y); 
    } 
    public Point ToTmsTile(Point t, Int32 zoom) { 
     return new Point(t.X, ((Int32)Math.Pow(2, zoom) - 1) - t.Y); 
    } 
} 
public struct Point { 
    public Point(Int32 x, Int32 y) 
     : this() { 
     X = x; 
     Y = y; 
    } 
    public Int32 X { get; set; } 
    public Int32 Y { get; set; } 
} 
public struct DPoint { 
    public DPoint(Double x, Double y) 
     : this() { 
     this.X = x; 
     this.Y = y; 
    } 
    public Double X { get; set; } 
    public Double Y { get; set; } 
    public static implicit operator DPoint(Point p) { 
     return new DPoint(p.X, p.Y); 
    } 
} 
public class GeoPoint { 
    public Double Latitude { get; set; } 
    public Double Longitude { get; set; } 
} 
public class Pair<T> { 
    public Pair() {} 
    public Pair(T min, T max) { 
     Min = min; 
     Max = max; 
    } 
    public T Min { get; set; } 
    public T Max { get; set; } 
} 

我有兩個問題。

1. 難道我口正確的代碼？ （我故意省略了，因爲我不使用它的一種方法，並添加一個我自己）

2. 這裏我有座標

 
WGS84 datum (longitude/latitude): 
-123.75 36.59788913307022 
-118.125 40.97989806962013 

Spherical Mercator (meters): 
-13775786.985667605 4383204.9499851465 
-13149614.849955441 5009377.085697312 

Pixels 
2560 6144 2816 6400 

Google 
x:10, y:24, z:6 

TMS 
x:10, y:39, z:6 

QuadTree 
023010

我應該如何鏈中的方法，使我從谷歌的獲得瓷磚座標（10,24,6）球形梅卡托米？

更新

找到適合我的第二個問題的答案是對我來說更重要。

Answer 1

有你的班上至少一個錯誤：

public Double Resolution(Int32 zoom) { 
     return InitialResolution/(2^zoom); // BAD CODE, USE Math.Pow, not^
    } 

如果你誤認爲是指數運算符二進制XOR運算符。

我已經重寫了代碼，使大多數功能靜態的，並且增加了一些相關的功能：

/// <summary> 
/// Conversion routines for Google, TMS, and Microsoft Quadtree tile representations, derived from 
/// http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/ 
/// </summary> 
public class WebMercator 
{ 
    private const int TileSize = 256; 
    private const int EarthRadius = 6378137; 
    private const double InitialResolution = 2 * Math.PI * EarthRadius/TileSize; 
    private const double OriginShift = 2 * Math.PI * EarthRadius/2; 

    //Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913 
    public static Point LatLonToMeters(double lat, double lon) 
    { 
     var p = new Point(); 
     p.X = lon * OriginShift/180; 
     p.Y = Math.Log(Math.Tan((90 + lat) * Math.PI/360))/(Math.PI/180); 
     p.Y = p.Y * OriginShift/180; 
     return p; 
    } 

    //Converts XY point from (Spherical) Web Mercator EPSG:3785 (unofficially EPSG:900913) to lat/lon in WGS84 Datum 
    public static Point MetersToLatLon(Point m) 
    { 
     var ll = new Point(); 
     ll.X = (m.X/OriginShift) * 180; 
     ll.Y = (m.Y/OriginShift) * 180; 
     ll.Y = 180/Math.PI * (2 * Math.Atan(Math.Exp(ll.Y * Math.PI/180)) - Math.PI/2); 
     return ll; 
    } 

    //Converts pixel coordinates in given zoom level of pyramid to EPSG:900913 
    public static Point PixelsToMeters(Point p, int zoom) 
    { 
     var res = Resolution(zoom); 
     var met = new Point(); 
     met.X = p.X * res - OriginShift; 
     met.Y = p.Y * res - OriginShift; 
     return met; 
    } 

    //Converts EPSG:900913 to pyramid pixel coordinates in given zoom level 
    public static Point MetersToPixels(Point m, int zoom) 
    { 
     var res = Resolution(zoom); 
     var pix = new Point(); 
     pix.X = (m.X + OriginShift)/res; 
     pix.Y = (m.Y + OriginShift)/res; 
     return pix; 
    } 

    //Returns a TMS (NOT Google!) tile covering region in given pixel coordinates 
    public static Tile PixelsToTile(Point p) 
    { 
     var t = new Tile(); 
     t.X = (int)Math.Ceiling(p.X/(double)TileSize) - 1; 
     t.Y = (int)Math.Ceiling(p.Y/(double)TileSize) - 1; 
     return t; 
    } 

    public static Point PixelsToRaster(Point p, int zoom) 
    { 
     var mapSize = TileSize << zoom; 
     return new Point(p.X, mapSize - p.Y); 
    } 

    //Returns tile for given mercator coordinates 
    public static Tile MetersToTile(Point m, int zoom) 
    { 
     var p = MetersToPixels(m, zoom); 
     return PixelsToTile(p); 
    } 

    //Returns bounds of the given tile in EPSG:900913 coordinates 
    public static Rect TileBounds(Tile t, int zoom) 
    { 
     var min = PixelsToMeters(new Point(t.X * TileSize, t.Y * TileSize), zoom); 
     var max = PixelsToMeters(new Point((t.X + 1) * TileSize, (t.Y + 1) * TileSize), zoom); 
     return new Rect(min, max); 
    } 

    //Returns bounds of the given tile in latutude/longitude using WGS84 datum 
    public static Rect TileLatLonBounds(Tile t, int zoom) 
    { 
     var bound = TileBounds(t, zoom); 
     var min = MetersToLatLon(new Point(bound.Left, bound.Top)); 
     var max = MetersToLatLon(new Point(bound.Right, bound.Bottom)); 
     return new Rect(min, max); 
    } 

    //Resolution (meters/pixel) for given zoom level (measured at Equator) 
    public static double Resolution(int zoom) 
    { 
     return InitialResolution/(Math.Pow(2, zoom)); 
    } 

    public static double ZoomForPixelSize(double pixelSize) 
    { 
     for (var i = 0; i < 30; i++) 
      if (pixelSize > Resolution(i)) 
       return i != 0 ? i - 1 : 0; 
     throw new InvalidOperationException(); 
    } 

    // Switch to Google Tile representation from TMS 
    public static Tile ToGoogleTile(Tile t, int zoom) 
    { 
     return new Tile(t.X, ((int)Math.Pow(2, zoom) - 1) - t.Y); 
    } 

    // Switch to TMS Tile representation from Google 
    public static Tile ToTmsTile(Tile t, int zoom) 
    { 
     return new Tile(t.X, ((int)Math.Pow(2, zoom) - 1) - t.Y); 
    } 

    //Converts TMS tile coordinates to Microsoft QuadTree 
    public static string QuadTree(int tx, int ty, int zoom) 
    { 
     var quadtree = ""; 

     ty = ((1 << zoom) - 1) - ty; 
     for (var i = zoom; i >= 1; i--) 
     { 
      var digit = 0; 

      var mask = 1 << (i - 1); 

      if ((tx & mask) != 0) 
       digit += 1; 

      if ((ty & mask) != 0) 
       digit += 2; 

      quadtree += digit; 
     } 

     return quadtree; 
    } 

    //Converts a quadtree to tile coordinates 
    public static Tile QuadTreeToTile(string quadtree, int zoom) 
    { 
     int tx= 0, ty = 0; 

     for (var i = zoom; i >= 1; i--) 
     { 
      var ch = quadtree[zoom - i]; 
      var mask = 1 << (i - 1); 

      var digit = ch - '0'; 

      if ((digit & 1) != 0) 
       tx += mask; 

      if ((digit & 2) != 0) 
       ty += mask; 
     } 

     ty = ((1 << zoom) - 1) - ty; 

     return new Tile(tx, ty); 
    } 

    //Converts a latitude and longitude to quadtree at the specified zoom level 
    public static string LatLonToQuadTree(Point latLon, int zoom) 
    { 
     Point m = LatLonToMeters(latLon.Y, latLon.X); 
     Tile t = MetersToTile(m, zoom); 
     return QuadTree(t.X, t.Y, zoom); 
    } 

    //Converts a quadtree location into a latitude/longitude bounding rectangle 
    public static Rect QuadTreeToLatLon(string quadtree) 
    { 
     int zoom = quadtree.Length; 
     Tile t = QuadTreeToTile(quadtree, zoom); 
     return TileLatLonBounds(t, zoom); 
    } 

    //Returns a list of all of the quadtree locations at a given zoom level within a latitude/longude box 
    public static List<string> GetQuadTreeList(int zoom, Point latLonMin, Point latLonMax) 
    { 
     if (latLonMax.Y< latLonMin.Y|| latLonMax.X< latLonMin.X) 
      return null; 

     Point mMin = LatLonToMeters(latLonMin.Y, latLonMin.X); 
     Tile tmin = MetersToTile(mMin, zoom); 
     Point mMax = LatLonToMeters(latLonMax.Y, latLonMax.X); 
     Tile tmax = MetersToTile(mMax, zoom); 

     var arr = new List<string>(); 

     for (var ty = tmin.Y; ty <= tmax.Y; ty++) 
     { 
      for (var tx = tmin.X; tx <= tmax.X; tx++) 
      { 
       var quadtree = QuadTree(tx, ty, zoom); 
       arr.Add(quadtree); 
      } 
     } 
     return arr; 
    } 
} 


/// <summary> 
/// Reference to a Tile X, Y index 
/// </summary> 
public class Tile 
{ 
    public Tile() { } 
    public Tile(int x, int y) 
    { 
     X = x; 
     Y = y; 
    } 

    public int X { get; set; } 
    public int Y { get; set; } 
} 

爲了解決第二個問題，請按照以下順序：

 int zoom = 6; 
     Tile googleTile = new Tile(10,24); 
     Tile tmsTile = GlobalMercator.ToTmsTile(googleTile, zoom); 
     Rect r3 = GlobalMercator.TileLatLonBounds(tmsTile, zoom); 
     var tl = GlobalMercator.LatLonToMeters(r3.Top, r3.Left); 
     var br = GlobalMercator.LatLonToMeters(r3.Bottom, r3.Right); 

     Debug.WriteLine("{0:0.000} {1:0.000}", tl.X, tl.Y); 
     Debug.WriteLine("{0:0.000} {1:0.000}", br.X, br.Y); 

     // -13775787.000 4383205.000 
     // -13149615.000 5009376.500 

Answer 2

從一個投影轉換座標到另一個最好的開源解決方案是Proj4原文爲C，但移植到衆多的編程語言。我試過並使用過的C＃端口是CodePlex上的DotSpatial Projections。基於示例很容易找到如何使用它。您需要知道的唯一事情就是您的案例的轉換參數。

Answer 3

一對於任何想要使用Oybek優秀代碼的讀者來說，指針的幾個：

你需要添加using System.Windows，但也Add a Reference th ËWindowsBase裝配，否則VS不會找到Point和Rect。

注意只是絕不能使用System.Drawing

這裏是一個新的功能，將變焦緯度/經度轉換成谷歌瓷磚：

public static Tile LatLonToGoogleTile(Point latLon, int zoom) 
    { 
     Point m = LatLonToMeters(latLon.Y, latLon.X); 
     Tile t = MetersToTile(m, zoom); 
     return ToGoogleTile(t, zoom); 
    }