国际上通用的是wgs84坐标系，而我国对于境内的坐标进行了加密，采用了gcj02坐标系，或者称为火星坐标系。亢孟军老师带的一门课《多媒体电子地图设计》要求我们从wgs84坐标系转换为gcj02坐标系，再反算出wgs84坐标系并进行前后wgs84坐标系误差分析。在这里简单介绍一下方法。
生成等间距矩阵点
首先简单介绍一段代码，使用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中做的，如下图所示，颜色越深代表误差越大，越浅则误差越小。
从图中可以看出中国的三级阶梯分布。