前文提供中国电子地图校正精确校正数据免费下载 提供了用于校正中国地图偏移的数据,下需给出如何利用这些数据来校正地图偏移的算法代码。
这个校正适用于 Google map China, Microsoft map china ,MapABC 等,因为这些地图构成方法是一样的。
有兴趣的可以参见 http://www.mapdigit.com/forum/vi ... &extra=page%3D1
具体原理这里不详细说,有兴趣的可以参见后续有关 Mapdigit 地图开发包教程或是 http://www.pstreets.com/gis.php 中开发指南。
关于校正精度的问题,0.01的数据相邻两个记录差最大在5个点左右,对应于18级每点代表0.5972米,也就是校正好误差在3米以内,而GPS接收器本身误差也在5米左右,比0.01再高的数据已无多大实际意义。5M以内对于精确定位及导航已足够,特别是对于导航时一般还会同时使用路径匹配算法,可以将位置精确校正路径上。
第一部分代码来自 http://msdn.microsoft.com/en-us/library/bb259689.aspx
第二部分代码为校正算法,算法很简单
view plaincopy to clipboardprint? using System; using System.Collections.Generic; using System.Linq; using System.Text; namespace MapDigit.GIS { public class GeoLatLng { public GeoLatLng(double latitude, double longitude) { this.latitude = latitude; this.longitude = longitude; } public double latitude; public double longitude; } public class GeoPoint { public GeoPoint(int x, int y) { this.x = x; this.y = y; } public int x; public int y; } public class OffsetInChina { //用于从GPS坐标转换为偏移后坐标 public static GeoLatLng fromEarthToMars(GeoLatLng earth) { GeoPoint ptOffset = getOffset(earth.latitude, earth.longitude); if (ptOffset.x != 0 || ptOffset.y != 0) { int pixelX, pixelY; TileSystem.LatLongToPixelXY(earth.latitude, earth.longitude, 18, out pixelX, out pixelY); GeoPoint pt = new GeoPoint(pixelX, pixelY); pt.x += ptOffset.x; pt.y += ptOffset.y; double latitude, longitude; TileSystem.PixelXYToLatLong(pt.x, pt.y, 18, out latitude, out longitude); return new GeoLatLng(latitude, longitude); } else { return new GeoLatLng(earth.latitude, earth.longitude); } } //用于将偏移后坐标转成真实的坐标 public static GeoLatLng fromMarToEarth(GeoLatLng mars) { GeoPoint ptOffset = getOffset(mars.latitude, mars.longitude); if (ptOffset.x != 0 || ptOffset.y != 0) { int pixelX, pixelY; TileSystem.LatLongToPixelXY(mars.latitude, mars.longitude, 18, out pixelX, out pixelY); GeoPoint pt = new GeoPoint(pixelX, pixelY); pt.x -= ptOffset.x; pt.y -= ptOffset.y; double latitude, longitude; TileSystem.PixelXYToLatLong(pt.x, pt.y, 18, out latitude, out longitude); return new GeoLatLng(latitude, longitude); } else { return new GeoLatLng(mars.latitude, mars.longitude); } } //这个函数用于将需要查询的经纬度转成最近的0.01分度值,无插值 //也可以自行实现插值 private static GeoPoint getQueryLocation(double latitude, double longitude) { int lat = (int)(latitude * 100); int lng = (int)(longitude * 100); double lat1 = ((int)(latitude * 1000 + 0.499999)) / 10.0; double lng1 = ((int)(longitude * 1000 + 0.499999)) / 10.0; for (double x = longitude; x < longitude + 1; x += 0.5) { for (double y = latitude; x < latitude + 1; y += 0.5) { if (x <= lng1 && lng1 < (x + 0.5) && lat1 >= y && lat1 < (y + 0.5)) { return new GeoPoint((int)(x + 0.5), (int)(y + 0.5)); } } } return new GeoPoint(lng, lat); } private static GeoPoint getOffset(double longitude, double latitude) { //这个函数用于返回查询结果,就是从校正数据中返回18级时x,y方偏移 //可以自行实现 return null; } } } using System; using System.Collections.Generic; using System.Linq; using System.Text;
namespace MapDigit.GIS { public class GeoLatLng {
public GeoLatLng(double latitude, double longitude) { this.latitude = latitude; this.longitude = longitude; } public double latitude; public double longitude; }
public class GeoPoint { public GeoPoint(int x, int y) { this.x = x; this.y = y; } public int x; public int y; }
public class OffsetInChina { //用于从GPS坐标转换为偏移后坐标 public static GeoLatLng fromEarthToMars(GeoLatLng earth) { GeoPoint ptOffset = getOffset(earth.latitude, earth.longitude); if (ptOffset.x != 0 || ptOffset.y != 0) { int pixelX, pixelY; TileSystem.LatLongToPixelXY(earth.latitude, earth.longitude, 18, out pixelX, out pixelY); GeoPoint pt = new GeoPoint(pixelX, pixelY); pt.x += ptOffset.x; pt.y += ptOffset.y; double latitude, longitude; TileSystem.PixelXYToLatLong(pt.x, pt.y, 18, out latitude, out longitude); return new GeoLatLng(latitude, longitude);
} else { return new GeoLatLng(earth.latitude, earth.longitude); }
}
//用于将偏移后坐标转成真实的坐标 public static GeoLatLng fromMarToEarth(GeoLatLng mars) { GeoPoint ptOffset = getOffset(mars.latitude, mars.longitude); if (ptOffset.x != 0 || ptOffset.y != 0) { int pixelX, pixelY; TileSystem.LatLongToPixelXY(mars.latitude, mars.longitude, 18, out pixelX, out pixelY); GeoPoint pt = new GeoPoint(pixelX, pixelY); pt.x -= ptOffset.x; pt.y -= ptOffset.y; double latitude, longitude; TileSystem.PixelXYToLatLong(pt.x, pt.y, 18, out latitude, out longitude); return new GeoLatLng(latitude, longitude);
} else { return new GeoLatLng(mars.latitude, mars.longitude); } }
//这个函数用于将需要查询的经纬度转成最近的0.01分度值,无插值 //也可以自行实现插值 private static GeoPoint getQueryLocation(double latitude, double longitude) { int lat = (int)(latitude * 100); int lng = (int)(longitude * 100); double lat1 = ((int)(latitude * 1000 + 0.499999)) / 10.0; double lng1 = ((int)(longitude * 1000 + 0.499999)) / 10.0; for (double x = longitude; x < longitude + 1; x += 0.5) { for (double y = latitude; x < latitude + 1; y += 0.5) { if (x <= lng1 && lng1 < (x + 0.5) && lat1 >= y && lat1 < (y + 0.5)) { return new GeoPoint((int)(x + 0.5), (int)(y + 0.5)); } } } return new GeoPoint(lng, lat); }
private static GeoPoint getOffset(double longitude, double latitude) { //这个函数用于返回查询结果,就是从校正数据中返回18级时x,y方偏移 //可以自行实现 return null; }
} }
其中需自已实现的是private static GeoPoint getOffset(double longitude, double latitude)
就是从下载的数据中,根据经纬度查出对应的X,Y向偏差。
附第一部分代码,
view plaincopy to clipboardprint? using System; using System.Text; namespace MapDigit.GIS { static class TileSystem { private const double EARTH_RADIUS = 6378137; private const double MIN_LATITUDE = -85.05112878; private const double MAX_LATITUDE = 85.05112878; private const double MIN_LONGITUDE = -180; private const double MAX_LONGITUDE = 180; /// <summary> /// Clips a number to the specified minimum and maximum values. /// </summary> /// <param name="n">The number to clip.</param> /// <param name="minValue">Minimum allowable value.</param> /// <param name="maxValue">Maximum allowable value.</param> /// <returns>The clipped value.</returns> private static double Clip(double n, double minValue, double maxValue) { return Math.Min(Math.Max(n, minValue), maxValue); } /// <summary> /// Determines the map width and height (in pixels) at a specified level /// of detail. /// </summary> /// <param name="levelOfDetail">Level of detail, from 1 (lowest detail) /// to 23 (highest detail).</param> /// <returns>The map width and height in pixels.</returns> public static uint MapSize(int levelOfDetail) { return (uint)256 << levelOfDetail; } /// <summary> /// Determines the ground resolution (in meters per pixel) at a specified /// latitude and level of detail. /// </summary> /// <param name="latitude">Latitude (in degrees) at which to measure the /// ground resolution.</param> /// <param name="levelOfDetail">Level of detail, from 1 (lowest detail) /// to 23 (highest detail).</param> /// <returns>The ground resolution, in meters per pixel.</returns> public static double GroundResolution(double latitude, int levelOfDetail) { latitude = Clip(latitude, MIN_LATITUDE, MAX_LATITUDE); return Math.Cos(latitude * Math.PI / 180) * 2 * Math.PI * EARTH_RADIUS / MapSize(levelOfDetail); } /// <summary> /// Determines the map scale at a specified latitude, level of detail, /// and screen resolution. /// </summary> /// <param name="latitude">Latitude (in degrees) at which to measure the /// map scale.</param> /// <param name="levelOfDetail">Level of detail, from 1 (lowest detail) /// to 23 (highest detail).</param> /// <param name="screenDpi">Resolution of the screen, in dots per inch.</param> /// <returns>The map scale, expressed as the denominator N of the ratio 1 : N.</returns> public static double MapScale(double latitude, int levelOfDetail, int screenDpi) { return GroundResolution(latitude, levelOfDetail) * screenDpi / 0.0254; } /// <summary> /// Converts a point from latitude/longitude WGS-84 coordinates (in degrees) /// into pixel XY coordinates at a specified level of detail. /// </summary> /// <param name="latitude">Latitude of the point, in degrees.</param> /// <param name="longitude">Longitude of the point, in degrees.</param> /// <param name="levelOfDetail">Level of detail, from 1 (lowest detail) /// to 23 (highest detail).</param> /// <param name="pixelX">Output parameter receiving the X coordinate in pixels.</param> /// <param name="pixelY">Output parameter receiving the Y coordinate in pixels.</param> public static void LatLongToPixelXY(double latitude, double longitude, int levelOfDetail, out int pixelX, out int pixelY) { latitude = Clip(latitude, MIN_LATITUDE, MAX_LATITUDE); longitude = Clip(longitude, MIN_LONGITUDE, MAX_LONGITUDE); double x = (longitude + 180) / 360; double sinLatitude = Math.Sin(latitude * Math.PI / 180); double y = 0.5 - Math.Log((1 + sinLatitude) / (1 - sinLatitude)) / (4 * Math.PI); uint mapSize = MapSize(levelOfDetail); pixelX = (int)Clip(x * mapSize + 0.5, 0, mapSize - 1); pixelY = (int)Clip(y * mapSize + 0.5, 0, mapSize - 1); } /// <summary> /// Converts a pixel from pixel XY coordinates at a specified level of detail /// into latitude/longitude WGS-84 coordinates (in degrees). /// </summary> /// <param name="pixelX">X coordinate of the point, in pixels.</param> /// <param name="pixelY">Y coordinates of the point, in pixels.</param> /// <param name="levelOfDetail">Level of detail, from 1 (lowest detail) /// to 23 (highest detail).</param> /// <param name="latitude">Output parameter receiving the latitude in degrees.</param> /// <param name="longitude">Output parameter receiving the longitude in degrees.</param> public static void PixelXYToLatLong(int pixelX, int pixelY, int levelOfDetail, out double latitude, out double longitude) { double mapSize = MapSize(levelOfDetail); double x = (Clip(pixelX, 0, mapSize - 1) / mapSize) - 0.5; double y = 0.5 - (Clip(pixelY, 0, mapSize - 1) / mapSize); latitude = 90 - 360 * Math.Atan(Math.Exp(-y * 2 * Math.PI)) / Math.PI; longitude = 360 * x; } /// <summary> /// Converts pixel XY coordinates into tile XY coordinates of the tile containing /// the specified pixel. /// </summary> /// <param name="pixelX">Pixel X coordinate.</param> /// <param name="pixelY">Pixel Y coordinate.</param> /// <param name="tileX">Output parameter receiving the tile X coordinate.</param> /// <param name="tileY">Output parameter receiving the tile Y coordinate.</param> public static void PixelXYToTileXY(int pixelX, int pixelY, out int tileX, out int tileY) { tileX = pixelX / 256; tileY = pixelY / 256; } /// <summary> /// Converts tile XY coordinates into pixel XY coordinates of the upper-left pixel /// of the specified tile. /// </summary> /// <param name="tileX">Tile X coordinate.</param> /// <param name="tileY">Tile Y coordinate.</param> /// <param name="pixelX">Output parameter receiving the pixel X coordinate.</param> /// <param name="pixelY">Output parameter receiving the pixel Y coordinate.</param> public static void TileXYToPixelXY(int tileX, int tileY, out int pixelX, out int pixelY) { pixelX = tileX * 256; pixelY = tileY * 256; } /// <summary> /// Converts tile XY coordinates into a QuadKey at a specified level of detail. /// </summary> /// <param name="tileX">Tile X coordinate.</param> /// <param name="tileY">Tile Y coordinate.</param> /// <param name="levelOfDetail">Level of detail, from 1 (lowest detail) /// to 23 (highest detail).</param> /// <returns>A string containing the QuadKey.</returns> public static string TileXYToQuadKey(int tileX, int tileY, int levelOfDetail) { StringBuilder quadKey = new StringBuilder(); for (int i = levelOfDetail; i > 0; i--) { char digit = '0'; int mask = 1 << (i - 1); if ((tileX & mask) != 0) { digit++; } if ((tileY & mask) != 0) { digit++; digit++; } quadKey.Append(digit); } return quadKey.ToString(); } /// <summary> /// Converts a QuadKey into tile XY coordinates. /// </summary> /// <param name="quadKey">QuadKey of the tile.</param> /// <param name="tileX">Output parameter receiving the tile X coordinate.</param> /// <param name="tileY">Output parameter receiving the tile Y coordinate.</param> /// <param name="levelOfDetail">Output parameter receiving the level of detail.</param> public static void QuadKeyToTileXY(string quadKey, out int tileX, out int tileY, out int levelOfDetail) { tileX = tileY = 0; levelOfDetail = quadKey.Length; for (int i = levelOfDetail; i > 0; i--) { int mask = 1 << (i - 1); switch (quadKey[levelOfDetail - i]) { case '0': break; case '1': tileX |= mask; break; case '2': tileY |= mask; break; case '3': tileX |= mask; tileY |= mask; break; default: throw new ArgumentException("Invalid QuadKey digit sequence."); } } } |
|