wgs84坐标系与gcj02坐标系转换误差分布图 | Mapping the Error in Transformation between WGS84 and GCJ02 Coordinations
2017-10-12 18:54
471 查看
国际上通用的是wgs84坐标系,而我国对于境内的坐标进行了加密,采用了gcj02坐标系,或者称为火星坐标系。亢孟军老师带的一门课《多媒体电子地图设计》要求我们从wgs84坐标系转换为gcj02坐标系,再反算出wgs84坐标系并进行前后wgs84坐标系误差分析。在这里简单介绍一下方法。
下面是源代码的类 LocationDivide,可以直接将这段代码拷贝进行使用。
具体的使用方法非常简单,如下:
即只要输入左下和右上角的经纬度坐标,再输入点之间的间隔,调用LocationDivide.compute_block()方法即可在划定的区域内,均匀生成等间距的采样点。需要使用时只需要调用LocationDivide.bounds变量即可,可以获得所有的点坐标,从左到右,从下到上依次排列。
这段代码非常实用,可以应用于多种场景,特别是需要生成格网时和需要均匀空间采样时可以使用。
该API可以将wgs84坐标系、百度坐标系等坐标系的点转换为gcj02坐标系中点的值。
特别注意的是,坐标点经纬度小数不超过6位,且一次最多传入40对坐标点。
wgs84和gcj02 相互转换JAVA代码,包括我的代码也主要使用了这位博主的代码:
http://www.cnblogs.com/94cool/p/4266907.html
windows phone上wgs84转换成gcj02的C#代码,可以通过这个反推出gcj02计算wgs84的方法。
https://on4wp7.codeplex.com/SourceControl/changeset/view/21483#353936
在这里还是附上我的代码,主要是使用了94cool博主的代码,修改为C#之后的代码,包括了两个类 Gps和PositionUtil:
调用方法直接使用PositionUtil.gcj_To_Gps84()等方法,传入相应的参数即可。
从图中可以看出中国的三级阶梯分布。
生成等间距矩阵点
首先简单介绍一段代码,使用python写的,作用是生成等间距矩阵点。下面是源代码的类 LocationDivide,可以直接将这段代码拷贝进行使用。
# Research region class LocationDivide(object): def __init__(self, bound, size): # minLat,minLon,maxLat,maxLon self.minLat = float(str(bound).split(',')[0]) self.minLon = float(str(bound).split(',')[1]) self.maxLat = float(str(bound).split(',')[2]) self.maxLon = float(str(bound).split(',')[3]) self.size = size # Seperate block into blocks def compute_block(self): if (self.maxLat - self.minLat) % self.size == 0: lat_count = (self.maxLat - self.minLat) / self.size else: lat_count = (self.maxLat - self.minLat) / self.size + 1 if (self.maxLon - self.minLon) % self.size == 0: lon_count = (self.maxLon - self.minLon) / self.size else: lon_count = (self.maxLon - self.minLon) / self.size + 1 self.bounds = [] lat_count = int(lat_count) lon_count = int(lon_count) try: for i in range(0, lat_count): for j in range(0, lon_count): # maxLat,minLon,minLat,maxLon minLat = self.minLat + i * self.size minLon = self.minLon + j * self.size maxLat = self.minLat + (i + 1) * self.size if maxLat > self.maxLat: maxLat = self.maxLat maxLon = self.minLon + (j + 1) * self.size if maxLon > self.maxLon: maxLon = self.maxLon # minLat,minLon,maxLat,maxLon # set decimal digit_number = 5 minLat = round(minLat, digit_number) minLon = round(minLon, digit_number) maxLat = round(maxLat, digit_number) maxLon = round(maxLon, digit_number) bound = "{0},{1},{2},{3}".format(minLat, minLon, maxLat, maxLon) self.bounds.append(bound) except Exception as e: with open("e:log.txt", 'a') as log: log.writelines(e) return self.bounds
具体的使用方法非常简单,如下:
# Set region bound and interval # minLat,minLon,maxLat,maxLon,interval region = "17,73,53,135" location = LocationDivide(region, 0.5) # Seperate grid into blocks location.compute_block()
即只要输入左下和右上角的经纬度坐标,再输入点之间的间隔,调用LocationDivide.compute_block()方法即可在划定的区域内,均匀生成等间距的采样点。需要使用时只需要调用LocationDivide.bounds变量即可,可以获得所有的点坐标,从左到右,从下到上依次排列。
这段代码非常实用,可以应用于多种场景,特别是需要生成格网时和需要均匀空间采样时可以使用。
wgs84坐标系转换为gcj02坐标系
在这里我们可以使用高德地图API进行数据的转换。具体可以参看其API该API可以将wgs84坐标系、百度坐标系等坐标系的点转换为gcj02坐标系中点的值。
特别注意的是,坐标点经纬度小数不超过6位,且一次最多传入40对坐标点。
wgs84坐标系与gcj02坐标系之间的相互转换
虽然我们可以使用高德地图API将wgs84数据转换为gcj02坐标系的数据,但如何从gcj02反解出wgs84坐标是一个问题。当然网上有很多相应的资料,这里仅仅简单列举两条我的参考文献:wgs84和gcj02 相互转换JAVA代码,包括我的代码也主要使用了这位博主的代码:
http://www.cnblogs.com/94cool/p/4266907.html
windows phone上wgs84转换成gcj02的C#代码,可以通过这个反推出gcj02计算wgs84的方法。
https://on4wp7.codeplex.com/SourceControl/changeset/view/21483#353936
在这里还是附上我的代码,主要是使用了94cool博主的代码,修改为C#之后的代码,包括了两个类 Gps和PositionUtil:
public class Gps { double latitude { set; get; } double longitude { set; get; } public Gps(double latitude, double lontitude) { this.latitude = latitude; this.longitude = lontitude; } public Gps(string gps) { this.latitude = Convert.ToDouble(gps.Split(',')[1]); this.longitude = Convert.ToDouble(gps.Split(',')[0]); } public double getWgLat() { return this.latitude; } public double getWgLon() { return this.longitude; } } public class PositionUtil { public static String BAIDU_LBS_TYPE = "bd09ll"; public static double pi = 3.1415926535897932384626; public static double a = 6378245.0; public static double ee = 0.00669342162296594323; /** * 84 to 火星坐标系 (GCJ-02) World Geodetic System ==> Mars Geodetic System * * @param lat * @param lon * @return */ public static Gps gps84_To_Gcj02(double lat, double lon) { if (outOfChina(lat, lon)) { return null; } double dLat = transformLat(lon - 105.0, lat - 35.0); double dLon = transformLon(lon - 105.0, lat - 35.0); double radLat = lat / 180.0 * pi; double magic = Math.Sin(radLat); magic = 1 - ee * magic * magic; double sqrtMagic = Math.Sqrt(magic); dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi); dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * pi); double mgLat = lat + dLat; double mgLon = lon + dLon; return new Gps(mgLat, mgLon); } /** * * 火星坐标系 (GCJ-02) to 84 * * @param lon * @param lat * @return * */ public static Gps gcj_To_Gps84(double lat, double lon) { Gps gps = transform(lat, lon); double lontitude = lon * 2 - gps.getWgLon(); double latitude = lat * 2 - gps.getWgLat(); return new Gps(latitude, lontitude); } /** * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法 将 GCJ-02 坐标转换成 BD-09 坐标 * * @param gg_lat * @param gg_lon */ public static Gps gcj02_To_Bd09(double gg_lat, double gg_lon) { double x = gg_lon, y = gg_lat; double z = Math.Sqrt(x * x + y * y) + 0.00002 * Math.Sin(y * pi); double theta = Math.Atan2(y, x) + 0.000003 * Math.Cos(x * pi); double bd_lon = z * Math.Cos(theta) + 0.0065; double bd_lat = z * Math.Sin(theta) + 0.006; return new Gps(bd_lat, bd_lon); } /** * * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法 * * 将 BD-09 坐标转换成GCJ-02 坐标 * * @param * bd_lat * @param bd_lon * @return */ public static Gps bd09_To_Gcj02(double bd_lat, double bd_lon) { double x = bd_lon - 0.0065, y = bd_lat - 0.006; double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * pi); double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * pi); double gg_lon = z * Math.Cos(theta); double gg_lat = z * Math.Sin(theta); return new Gps(gg_lat, gg_lon); } /** * (BD-09)-->84 * @param bd_lat * @param bd_lon * @return */ public static Gps bd09_To_Gps84(double bd_lat, double bd_lon) { Gps gcj02 = PositionUtil.bd09_To_Gcj02(bd_lat, bd_lon); Gps map84 = PositionUtil.gcj_To_Gps84(gcj02.getWgLat(), gcj02.getWgLon()); return map84; } public static Boolean outOfChina(double lat, double lon) { if (lon < 72.004 || lon > 137.8347) return true; if (lat < 0.8293 || lat > 55.8271) return true; return false; } public static Gps transform(double lat, double lon) { if (outOfChina(lat, lon)) { return new Gps(lat, lon); } double dLat = transformLat(lon - 105.0, lat - 35.0); double dLon = transformLon(lon - 105.0, lat - 35.0); double radLat = lat / 180.0 * pi; double magic = Math.Sin(radLat); magic = 1 - ee * magic * magic; double sqrtMagic = Math.Sqrt(magic); dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi); dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * pi); double mgLat = lat + dLat; double mgLon = lon + dLon; return new Gps(mgLat, mgLon); } public static 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 * Math.Sqrt(Math.Abs(x)); ret += (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0; ret += (20.0 * Math.Sin(y * pi) + 40.0 * Math.Sin(y / 3.0 * pi)) * 2.0 / 3.0; ret += (160.0 * Math.Sin(y / 12.0 * pi) + 320 * Math.Sin(y * pi / 30.0)) * 2.0 / 3.0; return ret; } public static double transformLon(double x, double y) { double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Math.Sqrt(Math.Abs(x)); ret += (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0; ret += (20.0 * Math.Sin(x * pi) + 40.0 * Math.Sin(x / 3.0 * pi)) * 2.0 / 3.0; ret += (150.0 * Math.Sin(x / 12.0 * pi) + 300.0 * Math.Sin(x / 30.0 * pi)) * 2.0 / 3.0; return ret; } }
调用方法直接使用PositionUtil.gcj_To_Gps84()等方法,传入相应的参数即可。
误差分布图
最后的误差分布图是在ArcGIS中做的,如下图所示,颜色越深代表误差越大,越浅则误差越小。从图中可以看出中国的三级阶梯分布。
相关文章推荐
- 百度坐标(BD09)、国测局坐标(火星坐标,GCJ02)、和WGS84坐标系之间的转换
- 百度地图与谷歌地图坐标转换,WGS84、GCJ02、BD09地图坐标系间的坐标转换及坐标距离计算
- WGS84、GCJ02、BD09各坐标系之间的转换算法
- WGS84 GCJ02和BD09坐标系相互转换代码
- 百度坐标(BD09)、国测局坐标(火星坐标,GCJ02)、和WGS84坐标系之间的转换(JS版代码)
- 一个提供了百度坐标(BD09)、国测局坐标(火星坐标,GCJ02)、和WGS84坐标系之间的转换的工具模块。
- 百度坐标(BD09)、国测局坐标(火星坐标,GCJ02)、和WGS84坐标系之间的转换
- 各地图坐标系转换(WGS84坐标系,GCJ02坐标系,BD09坐标系)
- 地球坐标系(WGS84),火星坐标系(GCJ02), 百度坐标系(BD09)坐标转换
- WGS84、GCJ02、BD09地图坐标系间的坐标转换及坐标距离计算
- [JS] 百度坐标(BD09)、国测局坐标(火星坐标,GCJ02)、和WGS84坐标系之间的转换(JS版代码)
- 百度坐标(BD09)、国测局坐标(火星坐标,GCJ02)、和WGS84坐标系之间的转换
- 地球坐标系(WGS84),火星坐标系(GCJ02), 百度坐标系(BD09)坐标转换
- 关于百度坐标系 (BD-09)与火星坐标系 (GCJ-02)以及WGS84坐标之间的互相转换
- WGS84坐标系和CGS2000国家坐标系统转换
- 北京54全国80及WGS84坐标系的相互转换
- Google投影坐标和WGS84坐标系的转换
- 如何将谷歌地图高清卫星影像坐标系转换成西安80坐标(WGS84)
- 关于GCJ02和WGS84坐标系的一点实验