Convertisseur de coordonnées géographiques

Question

Je voulais une classe qui puisse convertir un système en un autre.

J'ai trouvé un code source en python et j'ai essayé de le porter en C#.

Ceci est la source 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 

Voici mon port 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; } 
} 

J'ai deux questions.

1. Ai-je porté le code correctement? (I a intentionnellement omis une méthode comme je ne l'utilise pas et a ajouté un ma propre)

2. Ici, j'ai coordonnées

 
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

Comment dois-je enchaîner les méthodes, afin que je reçois de Google tuiles coordonne (10, 24, 6) les compteurs mercator sphériques?

Mise à jour

Trouver réponse à ma deuxième question est plus important pour moi.

Answer 1

Il y a au moins un bug dans votre classe:

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

Si vous avez trompé l'opérateur XOR binaire pour l'opérateur d'exposant.

J'ai réécrit le code, a fait la plupart des fonctions statiques, et a ajouté quelques fonctions plus pertinentes:

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

Pour résoudre votre deuxième question, utilisez la séquence suivante:

 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

La meilleure solution opensource pour convertir les coordonnées d'une projection à une autre est Proj4 écrit à l'origine en c mais porté vers de nombreux langages de programmation. Le port à C# que j'ai essayé et utilisé est DotSpatial Projections trouvé sur CodePlex. Il est facile de savoir comment l'utiliser en fonction des exemples. La seule chose que vous devez savoir sont les paramètres de conversion pour votre cas.

Answer 3

A couple de pointeurs pour toute lecture qui veut utiliser Oybek excellent code:

Vous devez ajouter using System.Windows, mais aussi Add a Reference à th e WindowsBase assemblage, sinon VS ne trouverez pas Point et Rect.

Notez que juste ne doit pas utiliserSystem.Drawing

Et voici une nouvelle fonction qui convertira Zoom lat/lng à une tuile 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); 
    }