地理坐标转换器

Question

我需要一个能够将一个系统转换为另一个系统的类。

我找到了一份Python源代码并尝试将其移植到C#中。

这是Python源代码。 从这里获取。

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基准（经度/纬度）：
-123.75 36.59788913307022
-118.125 40.97989806962013
球形墨卡托（米）：
-13775786.985667605 4383204.9499851465
-13149614.849955441 5009377.085697312
像素：
2560 6144 2816 6400
Google：
x:10, y:24, z:6
TMS：
x:10, y:39, z:6
QuadTree：
023010

我应该如何链接这些方法，以便从谷歌的瓦片坐标（10、24、6）获取球形墨卡托米数？

更新

对于我的第二个问题，找到答案对我更重要。

- Oybek

---

只需运行并测试它。如果您获得所需的输出，那么您可能已经成功了。 - Rohit Jain

我对GIS并不完全精通，因此我正在寻找专家的帮助。 - Oybek

自从这篇文章发布以来，链接页面底部已经有了源Python代码的完整C#端口。 - alofgran

Answer 1

你的类中至少有一个 bug：

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

你把二进制异或运算符误认为指数运算符了。

我已经重写了代码，使大部分函数静态化，并增加了一些相关的函数。

/// <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

- Jay Borseth

---

我认为你的 MetersToLatLon 函数有误。你把 X 和 Y 值搞反了。 - Tri Q Tran

"GlobalMercator"是什么？它是“Web Mercator”[https://en.wikipedia.org/wiki/Web_Mercator]还是“Sphere Mercator”[https://en.wikipedia.org/wiki/Mercator_projection]？ - Kind Contributor

根据上面提到的Python原始来源，GlobalMercator是EPSG 3785（后来成为EPSG 3857）。Google地图是基于WGS84（EPSG 3857）的投影坐标系统。此外，阅读此链接以获取更多信息：https://gis.stackexchange.com/questions/48949/epsg-3857-or-4326-for-googlemaps-openstreetmap-and-leaflet - Robert Smith

尽管代码无法构建，Point类未定义，并且Point和Tile类似乎是多余的，但这让我开始了将经纬度坐标转换为X/Y坐标，以便在请求地图瓦片时使用四叉树键。 - Justin

尽管代码无法构建，Point类未定义，而Point和Tile类似乎是多余的，但这让我开始将经纬度坐标转换为X/Y坐标，以便在请求地图瓦片的四叉树键中使用。 - Justin

Answer 2

将坐标从一个投影转换为另一个的最佳开源解决方案是Proj4，它最初是用c编写的，但被移植到了许多编程语言中。我尝试并使用的c#移植版本是DotSpatial Projections，可以在CodePlex上找到。它很容易根据示例了解如何使用。唯一需要知道的是您的情况下的转换参数。

- Nikola Davidovic

---

链接已失效：http://dotspatial.codeplex.com/wikipage?title=DotSpatial.Projections - Gino Mempin

这个链接似乎有效：https://github.com/DotSpatial/DotSpatial，但是它指向更深层次的许多失效链接... - joanis

Answer 3

CoordinateSharp 可以在NuGet上获取。它非常轻巧，可以实现坐标转换。它不是为地图设计的，而是针对纯粹的转换（和基于位置的天文信息），如果这是一个因素。

示例

 Coordinate c = new Coordinate(myLat,myLong);
 c.UTM; //Outputs UTM string
 c.UTM.Easting //UTM easting property

- Tronald

Answer 4

对于想要使用Oybek优秀代码的读者，这里有几个提示：

你需要添加using System.Windows，但也需要添加引用到WindowsBase程序集，否则VS将无法找到Point和Rect。

请注意绝不能使用System.Drawing

以下是一个新函数，可以将缩放的经纬度转换为Google瓦片：

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

- NickPR