Geographic coordinates converter

Question

I wanted a class that can convert one system to another.

I've found a source code in python and tried to port it into C#.

This is the python source. 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

Here is my C# port.

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; }
}

I have two questions.

1. Did I port the code correctly? (I intentionally omitted one method as I don't use it and added one my own)

2. Here I have coordinates

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

How should I chain the methods so that I get from Google's tile coordinates (10, 24, 6) the spherical mercator meters?

Update

Finding answer for my second question is more important for me.

Answer 1

There's at least one bug in your class:

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

Where you've mistaken the binary XOR operator for the exponent operator.

I've rewritten the code, made most functions static, and added a few more relevant functions:

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

To solve your second question, use the following sequence:

        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

The best opensource solution for converting coordinates from one projection to another is Proj4 originally written in c but ported to numerous programming languages. The port to c# that I have tried and used is DotSpatial Projections found on CodePlex. It is easy to find out how to use it based on the examples. The only thing you need to know are conversion parameters for your case.