中国地图偏移算法.doc

WGS-84 坐标系，是 GPS 系统所采用的坐标系。一切正常工作的 GPS 或 GPS 芯片所返回的坐标值都是这个坐标系下的数值。Google 地图采用的卫星图也是按照这个坐标系摆放的。GCJ-02 坐标系，是我天朝政府搞出来的加密坐标系，也常常被称为 “火星坐标系”。包括（但可能不限于）高德地图在内的国内地图服务商采用它来绘制地图。Apple、Google 等国外公司在其道路地图中使用的也是高德的数据。BD-09 坐标系则是百度地图专用的坐标系。在开发一些 LBS 应用时，如果不加处理，很容易出现几种形式的地图之间出现偏移的情况。因此在这几个坐标系之间进行转换非常重要。以下代码就是网络上泄露出的从 WGS-84 转 GCJ-02 的算法，以及 GCJ-02 与 BD-09 的互相转换算法。另外，代码中的从 GCJ-02 到 WGS-84 的逆变换是我用迭代法实现的。实验证明，由于原变换的局部线性性非常好，只要迭代两次（调用原变换三次），经纬度就可以达到大约 1e-8 度的精度（相当于毫米级精度）。所以运行速度不会太慢，对于一般应用足矣。全部代码也可以在 GitHub 下载：/fengzee-me/ChinaMapShift。代码和其编译出的 libChinaMapShift.a 库文件直接适用于 iOS 环境下，需要在其他平台下运行时，可以自行编译（源代码完全是原生 C 代码，只需要 math.h 库）。12345678910111213141516171819202122/ ChinaMapShift.h/ ChinaMapShift/ Most code created by someone anonymous./ transformFromGCJToWGS() added by Fengzee ()./ #ifndef ChinaMapShift_ChinaMapShift_h#define ChinaMapShift_ChinaMapShift_h typedef struct double lng; double lat; Location; Location transformFromWGSToGCJ(Location wgLoc);Location transformFromGCJToWGS(Location gcLoc);Location bd_encrypt(Location gcLoc);Location bd_decrypt(Location bdLoc); #endif中国地图偏移算法源文件C123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129/ ChinaMapShift.mm/ ChinaMapShift/ Most code created by someone anonymous./ transformFromGCJToWGS() added by Fengzee ()./ #include ChinaMapShift.h#include inline LocationLocationMake(double lng, double lat) Location loc; loc.lng = lng, loc.lat = lat; return loc; / Transform WGS-84 to GCJ-02 (Chinese encrypted coordination system)/ const double pi = 3.14159265358979324; / Krasovsky 1940/ a = 6378245.0, 1/f = 298.3/ b = a * (1 - f)/ ee = (a2 - b2) / a2;const double a = 6378245.0;const double ee = 0.00669342162296594323; int outOfChina(double lat, double lon) if (lon 137.8347) return 1; if (lat 55.8271) return 1; return 0; double transformLat(double x, double y) double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(x 0 ? x:-x); ret += (20.0 * sin(6.0 * x * pi) + 20.0 *sin(2.0 * x * pi) * 2.0 / 3.0; ret += (20.0 * sin(y * pi) + 40.0 * sin(y / 3.0 * pi) * 2.0 / 3.0; ret += (160.0 * sin(y / 12.0 * pi) + 320 * sin(y * pi / 30.0) * 2.0 / 3.0; return ret; double transformLon(double x, double y) double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(x 0 ? x:-x); ret += (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi) * 2.0 / 3.0; ret += (20.0 * sin(x * pi) + 40.0 * sin(x / 3.0 * pi) * 2.0 / 3.0; ret += (150.0 * sin(x / 12.0 * pi) + 300.0 * sin(x / 30.0 * pi) * 2.0 / 3.0; return ret; Location transformFromWGSToGCJ(Location wgLoc) Location mgLoc; if (outOfChina(wgLoc.lat, wgLoc.lng) mgLoc = wgLoc; return mgLoc; double dLat = transformLat(wgLoc.lng - 105.0, wgLoc.lat - 35.0); double dLon = transformLon(wgLoc.lng - 105.0, wgLoc.lat - 35.0); double radLat = wgLoc.lat / 180.0 * pi; double magic = sin(radLat); magic = 1 - ee * magic * magic; double sqrtMagic = sqrt(magic); dLat = (dLat * 180.0) / (a * (1 - ee) / (magic * sqrtMagic) * pi); dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * pi); mgLoc.lat = wgLoc.lat + dLat; mgLoc.lng = wgLoc.lng + dLon; return mgLoc; / Transform GCJ-02 to WGS-84/ Reverse of transformFromWGSToGC() by iteration./ Created by Fengzee () on 13-11-9./Location transformFromGCJToWGS(Location gcLoc) Location wgLoc = gcLoc; Location currGcLoc, dLoc; while (1) currGcLoc = transformFromWGSToGCJ(wgLoc); dLoc.lat = gcLoc.lat - currGcLoc.lat; dLoc.lng = gcLoc.lng - currGcLoc.lng; if (fabs(dLoc.lat) 1e-7 & fabs(dLoc.lng) 1e-7) / 1e-7 centimeter level accuracy / Result of experiment: / Most of the time 2 iterations would be enough for an 1e-8 accuracy (milimeter level). / return wgLoc; wgLoc.lat += dLoc.lat; wgLoc.lng += dLoc.lng; return wgLoc; / Transform GCJ-02 to BD-09/const double x_pi = 3.14159265358979324 * 3000.0 / 180.0;Location bd_encrypt(Location gcLoc) double x = gcLoc.lng, y = gcLoc.lat; double z = sqrt(x * x + y * y) + 0.00002 * sin(y * x_pi); double theta = atan2(y, x) + 0.000003 * cos(x * x_pi); return LocationMake(z * cos(theta) + 0.0065, z * sin(theta) + 0.006); / Transform BD-09 to GCJ-02/Location bd_decrypt(Locat