Geographic coordinate converter

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

I found the source code in python and tried porting it to C #.

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

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(); pX = lon * OriginShift / 180; pY = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180); pY = pY * OriginShift / 180; return p; } public GeoPoint MetersToLatLon(DPoint m) { var ll = new GeoPoint(); ll.Longitude = (mX / OriginShift) * 180; ll.Latitude = (mY / 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 = pX * res - OriginShift; met.Y = pY * res - OriginShift; return met; } public DPoint MetersToPixels(DPoint m, Int32 zoom) { var res = Resolution(zoom); var pix = new DPoint(); pix.X = (mX + OriginShift) / res; pix.Y = (mY + OriginShift) / res; return pix; } public Point PixelsToTile(DPoint p) { var t = new Point(); tX = (Int32)Math.Ceiling(pX / (Double)TileSize) - 1; tY = (Int32)Math.Ceiling(pY / (Double)TileSize) - 1; return t; } public Point PixelsToRaster(Point p, Int32 zoom) { var mapSize = TileSize << zoom; return new Point(pX, mapSize - pY); } 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(tX * TileSize, tY * TileSize), zoom); var max = PixelsToMeters(new DPoint((tX + 1) * TileSize, (tY + 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(tX, ((Int32)Math.Pow(2, zoom) - 1) - tY); } public Point ToTmsTile(Point t, Int32 zoom) { return new Point(tX, ((Int32)Math.Pow(2, zoom) - 1) - tY); } } 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(pX, pY); } } 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.

  • Am I loading the code correctly? (I intentionally skipped one method since I don't use it, and added one of my own)

  • Here I have the 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 can I connect methods to get spherical mercantary counters from the coordinates of Google tiles (10, 24, 6)?

Update

Finding the answer to my second question is more important to me.

+4
source share
3 answers

There is at least one error in your class:

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

If you mistakenly accept the binary XOR operator for the exponent operator.

I rewrote the code, made most of the functions static and added some more important 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(); pX = lon * OriginShift / 180; pY = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180); pY = pY * 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 = (mX / OriginShift) * 180; ll.Y = (mY / 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 = pX * res - OriginShift; met.Y = pY * 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 = (mX + OriginShift) / res; pix.Y = (mY + 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(); tX = (int)Math.Ceiling(pX / (double)TileSize) - 1; tY = (int)Math.Ceiling(pY / (double)TileSize) - 1; return t; } public static Point PixelsToRaster(Point p, int zoom) { var mapSize = TileSize << zoom; return new Point(pX, mapSize - pY); } //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(tX * TileSize, tY * TileSize), zoom); var max = PixelsToMeters(new Point((tX + 1) * TileSize, (tY + 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(tX, ((int)Math.Pow(2, zoom) - 1) - tY); } // Switch to TMS Tile representation from Google public static Tile ToTmsTile(Tile t, int zoom) { return new Tile(tX, ((int)Math.Pow(2, zoom) - 1) - tY); } //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(tX, tY, 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 the 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 
+7
source

The best open source solution for converting coordinates from one projector to another is Proj4, originally written in c, but ported to many programming languages. The port in C # that I tried and used, DotSpatial Projections is found on CodePlex. It’s easy to learn how to use it based on examples. The only thing you need to know is the conversion options for your case.

+3
source

A few pointers for those who read that Oybek wants to use excellent code:

You need to add using System.Windows , as well as Add a Reference to the WindowsBase assembly, otherwise VS will not find Point and Rect .

Please note that just should not use System.Drawing

And here is a new feature that will convert Zoom lat / lng to Google Tile:

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

Source: https://habr.com/ru/post/1439818/


All Articles