c# - 地理座標コンバータ

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

2つの質問があります。

1. コードを正しく移植しましたか？（使用しないため、意図的に1つのメソッドを省略し、独自のメソッドを追加しました）

2. ここに座標があります

WGS84データム（経度/緯度）：
-123.75 36.59788913307022
-118.125 40.97989806962013

球形メルカトル図法（メートル）：
-13775786.985667605 4383204.9499851465
-13149614.849955441 5009377.085697312

ピクセル
2560 6144 2816 6400

グーグル
x：10、y：24、z：6

TMS
x：10、y：39、z：6

QuadTree
023010

Googleのタイル座標（10、24、6）から球形メルカトル図法を取得するには、メソッドをどのように連鎖させる必要がありますか？

アップデート

2番目の質問に対する答えを見つけることは私にとってより重要です。

Answer 1

クラスに少なくとも 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; }
}

2 番目の質問を解決するには、次のシーケンスを使用します。

        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

座標をある投影法から別の投影法に変換するための最良のオープンソース ソリューションは、もともと c で記述された Proj4 ですが、多数のプログラミング言語に移植されています。私が試して使用した c# への移植は、CodePlex にあるDotSpatial Projectionsです。例に基づいて使用方法を見つけるのは簡単です。知っておく必要があるのは、ケースの変換パラメーターだけです。

Answer 3

CoordinateSharpは NuGet で利用できます。軽量で、座標変換が本当にできます。マッピング用には設計されていませんが、それが要因である場合は、変換 (および位置に基づく天体情報) を制限します。

例

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

Answer 4

Oybek の優れたコードを使用したい読者のためのいくつかの指針:

追加する必要がありますが、アセンブリにusing System.Windowsも追加する必要があります。そうしないと、VS はと を見つけられません。Add a ReferenceWindowsBasePointRect

ただ使用してはならないことに注意してくださいSystem.Drawing

Zoom lat/lng を 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);
    }