• 1、地理坐标下切片系统的计算 地理坐标下切片系统的计算，主要适用于google地球中切片系统，以及目标底图参考系统为EPSG:4326的情况。 public class GlobalGeodetic { private int tileSize; private ...
1、地理坐标下切片系统的计算

public class GlobalGeodetic {
private int tileSize;
private double resFact;

public GlobalGeodetic(String tmscompatible, int tileSize) {
this.tileSize = tileSize;
if (tmscompatible != null && tmscompatible.length() > 0) {
// Defaults the resolution factor to 0.703125 (2 tiles @ level 0)
this.resFact = 180.0D / (double)this.tileSize;
} else {
//Defaults the resolution factor to 1.40625 (1 tile @ level 0)
this.resFact = 360.0D / (double)this.tileSize;
}

}

public double[] lonlatToPixels(double lon, double lat, int zoom) {
double res = this.resFact / Math.pow(2.0D, (double)zoom);
return new double[]{(180.0D + lon) / res, (90.0D + lat) / res};
}

public int[] pixelsToTile(double px, double py) {
int tx = (int)(Math.ceil(px / (double)this.tileSize) - 1.0D);
int ty = (int)(Math.ceil(py / (double)this.tileSize) - 1.0D);
return new int[]{tx, ty};
}

public int[] lonlatToTile(double lon, double lat, int zoom) {
double[] pxpy = this.lonlatToPixels(lon, lat, zoom);
return this.pixelsToTile(pxpy[0], pxpy[1]);
}

public double resolution(int zoom) {
return this.resFact / Math.pow(2.0D, (double)zoom);
}

public int zoomForPixelSize(double pixelSize) {
for(int i = 0; i < 32; ++i) {
if (pixelSize > this.resolution(i)) {
if (i != 0) {
return i - 1;
}

return 0;
}
}

return 0;
}

public double[] tileBounds(int tx, int ty, int zoom) {
double res = this.resFact / Math.pow(2.0D, (double)zoom);
return new double[]{(double)(tx * this.tileSize) * res - 180.0D, (double)(ty * this.tileSize) * res - 90.0D, (double)((tx + 1) * this.tileSize) * res - 180.0D, (double)((ty + 1) * this.tileSize) * res - 90.0D};
}

public double[] tileLatLonBounds(int tx, int ty, int zoom) {
double[] b = this.tileBounds(tx, ty, zoom);
return new double[]{b[1], b[0], b[3], b[2]};
}
}

2、投影坐标系下计算类

目前绝大多数切片系统内部实现都是采用这种方式实现。

public class GlobalMercator {
private int tileSize;
private double initialResolution;
private double originSh;
//6378137 为地球半径
public GlobalMercator(int tileSize) {
this.tileSize = tileSize;
this.initialResolution = 2 * Math.PI * 6378137 / this.tileSize;
this.originShift = 2 * Math.PI * 6378137 / 2.0;
}

public double[] latLonToMeters(double lat, double lon) {
double mx = lon * this.originShift / 180.0;
double my = Math.log(Math.tan((90 + lat) * Math.PI / 360.0)) / (Math.PI / 180.0);
my = my * this.originShift / 180.0;
return new double[]{mx, my};
}

public double[] metersToLatLon(double mx, double my) {
double lon = (mx / this.originShift) * 180.0;
double lat = (my / this.originShift) * 180.0;
lat = 180 / Math.PI * (2 * Math.atan(Math.exp(lat * Math.PI / 180.0)) - Math.PI / 2.0);
return new double[]{lat, lon};
}

public double[] pixelsToMeters(int px, int py, int zoom) {
double res = this.resolution(zoom);
double mx = px * res - this.originShift;
double my = py * res - this.originShift;
return new double[]{mx, my};
}

public double[] metersToPixels(double mx, double my, int zoom) {
double res = this.resolution(zoom);
double px = (mx + this.originShift) / res;
double py = (my + this.originShift) / res;
return new double[]{px, py};
}

public int[] pixelsToTile(double px, double py) {
int tx = (int) (Math.ceil(px / (float) (this.tileSize)) - 1);
int ty = (int) (Math.ceil(py / (float) (this.tileSize)) - 1);
return new int[]{tx, ty};
}

public double[] pixelsToRaster(double px, double py, int zoom) {
double mapSize = this.tileSize << zoom;
return new double[]{px, mapSize - py};
}

public int[] metersToTile(double mx, double my, int zoom) {
double[] coordinate = this.metersToPixels(mx, my, zoom);
return this.pixelsToTile(coordinate[0], coordinate[1]);
}

public double[] tileBounds(int tx, int ty, int zoom) {
double[] minxy = pixelsToMeters(tx * this.tileSize, ty * this.tileSize, zoom);
double[] maxxy = pixelsToMeters((tx + 1) * this.tileSize, (ty + 1) * this.tileSize, zoom);
return new double[]{minxy[0], minxy[1], maxxy[0], maxxy[1]};
}

public double[] tileLatLonBounds(int tx, int ty, int zoom) {
double[] bounds = this.tileBounds(tx, ty, zoom);
double[] minLatLon = this.metersToLatLon(bounds[0], bounds[1]);
double[] maxLatlon = this.metersToLatLon(bounds[2], bounds[3]);
return new double[]{minLatLon[0], minLatLon[1], maxLatlon[0], maxLatlon[1]};
}

public double resolution(int zoom) {
return this.initialResolution / (Math.pow(2, zoom));
}

public int zoomForPixelSize(double pixelSize) {
for (int i = 0; i < 32; i++) {
if (pixelSize > this.resolution(i)) {
if (i != 0) {
return i - 1;
} else return 0;
}
}
return 0;
}
public int[] googleTile(int tx, int ty, int zoom) {
return new int[]{tx, ((int) (Math.pow(2, zoom)) - 1) - ty};
}

public String quadTree(int tx, int ty, int zoom) {
ty = (int) (Math.pow(2, zoom) - 1 - ty);
for (int i = zoom; i > 0; i--) {
int digit = 0;
int mask = 1 << (i - 1);
if ((tx & mask) != 0) {
digit += 1;
}
if ((ty & mask) != 0) {
digit += 2;
}
}
}
//XYZ转换为geohash编码
public String tileXYToQuadKey(int tileX, int tileY, int levelOfDetail) {
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++;
}
}
}
//geohash编码转换为XYZ
int tileX = 0, tileY = 0;
for (int i = levelOfDetail; i > 0; i--) {
int mask = 1 << (i - 1);
case '0':
break;
case '1':
break;
case '2':
break;
case '3':
break;
default:
throw new Exception("Invalid QuadKey digit sequence.");
}
}
return new int[]{tileX, tileY, levelOfDetail};
}
}

3、使用例子

根据XYZ计算经纬度范围

//地理坐标（EPSG:4326）下计算方式
double[] bbox = new GlobalGeodetic("", 256).tileLatLonBounds(x, y, z);
//投影坐标（EPSG:3857）下的计算方式
double[] bboxs = new GlobalMercator(256).tileLatLonBounds(x, y, z);


根据经纬度范围计算XYZ

//地理坐标（EPSG:4326）下计算方式
……待补充

//投影坐标（EPSG:3857）下的计算方式
Envelope envelope = new Envelope(xmin, xmax, ymin, ymax);
GlobalMercator mercator = new GlobalMercator(256);
double[] min = mercator.latLonToMeters(envelope.getMinY(), envelope.getMinX());
double[] max = mercator.latLonToMeters(envelope.getMaxY(), envelope.getMaxX());

//#region 计算
for (int tz = tmaxz; tz > tminz - 1; tz--) {
int[] tminxy = mercator.metersToTile(min[0], min[1], tz);
int[] tmaxxy = mercator.metersToTile(max[0], max[1], tz);
tminxy = new int[]{Math.max(0, tminxy[0]), Math.max(0, tminxy[1])};
tmaxxy = new int[]{(int) Math.min(Math.pow(2, tz) - 1, tmaxxy[0]), (int) Math.min(Math.pow(2, tz) - 1, tmaxxy[1])};
for (int tx = tminxy[0]; tx < tmaxxy[0] + 1; tx++) {
for (int ty = tmaxxy[1]; ty > tminxy[1] - 1; ty--) {
//z,x,y坐标
}
}
}
//endregion



展开全文
• ## 地理坐标系与投影坐标系的区别

万次阅读 多人点赞 2018-08-17 22:57:17
平时开展GIS开发、研究、应用工作，总会接触到坐标系，也会遇到坐标转换的问题，如地理坐标系、投影坐标系等。 地理坐标系是球面坐标，参考平面是椭球面，坐标单位是经纬度； 投影坐标系是平面坐标系，参考平面...
1.基本概念

平时开展GIS开发、研究、应用工作，总会接触到坐标系，也会遇到坐标转换的问题，如地理坐标系、投影坐标系等。

地理坐标系是球面坐标，参考平面是椭球面，坐标单位是经纬度；

投影坐标系是平面坐标系，参考平面是水平面，坐标单位是米、千米等。

地理坐标系转换到投影坐标系的过程理解为投影，即将不规则的地球曲面转换为平面。

在当前的信息化的技术条件下，直接使用地理坐标系是不是更加真实准确，像谷歌地球；投影毕竟存在各种变形。

地理坐标系的WKID介绍：Geographic Coordinate Systems

投影坐标系的WKID介绍：Projected Coordinate Systems

EPSG:European Petroleum Survey Group，欧洲石油调查组织，

该组织负责专门维护地球上所有的测量坐标系统，并且给每组坐标系统都赋予了一个编号和一组描述（WKT），

比如大家常用的WGS84坐标系编号就是EPSG:4326，再比如互联网地图（谷歌、高德等）常用的伪墨卡托投影编号就是EPSG:3857。

可以理解成EPSG给大家维护了无数把尺子，并且给每把尺子搞了个编号，还标明了这把尺子适合什么条件下用。

2. 地理坐标系

2.1 地球的三级逼近

2.1.1 大地水准面

地球的自然表面不是平整的，需要想办法用数学公式描述地球表面，只能设想一个近似的数学面。

大地水准面是地球表面的第一级逼近。假设当海水处于完全静止的平衡状态时，从海平面延伸到所有大陆下部，而与地球重力方向处处正交的一个连续、闭合的曲面，这就是大地水准面。

地球椭球体是地球表面的第二级逼近。大地水准面可以近似成一个规则成椭球体，但并不是完全规则，其形状接近一个扁率极小的椭圆绕短轴旋转所形成的规则椭球体，这个椭球体称为地球椭球体。

地球椭球体的基本参数：

长半轴（赤道半径）    a
短半轴（极半径）      b
椭球体的扁率         à=(a-b)/a
第一偏心率           è=（a2-b2）/a2
第二偏心率           é=(a2-b2)/ b2


常见的椭球体的参数：

	克拉索夫斯基椭球	   1975 GRS椭球体	   WGS-84椭球体
a	6 378 245.000 m	   6 378 140.000 m	   6 378 137.000 m
b	6 356 863.019 m	   6 356 755.288 m	   6 356 752.314 m
à	   1/298.3	         1/298.257	        1/298.257 224
è	0.006 693 422	   0.006 694 385	   0.006 694 380
é	0.006 738 525	   0.006 739 502	   0.006 739 497


大地基准面是地球表面的第三极逼近。

椭球体是对地球的抽象，不能与地球表面完全重合，在设置参考椭球体的时候必然会出现有的地方贴近的好（参考椭球体与地球表面位置接近），有地地方贴近的不好的问题，因此这里还需要一个大地基准面来控制参考椭球和地球的相对位置。有以下两类基准面：

地心基准面：由卫星数据得到，使用地球的质心作为原点，使用最广泛的是 WGS 1984。

区域基准面：特定区域内与地球表面吻合，大地原点是参考椭球与大地水准面相切的点，例如Beijing-54、Xian-80。称谓的Beijing-54、Xian-80坐标系实际上指的是我国的两个大地基准面。

地心大地坐标系：指经过定位与定向后，地球椭球的中心与地球质心重合。如CGCS2000、WGS84。

参心大地坐标系：指经过定位与定向后，地球椭球的中心不与地球质心重合而是接近地球质心。区域性大地坐标系是我国基本测图和常规大地测量的基础。如Beijing-54、Xian-80。

2.2 地理坐标

地理坐标，就是用经线（子午线）、纬线、经度、纬度表示地面点位的球面坐标。

一般地理坐标可分为三种，天文经纬度，大地经纬度，地心经纬度。通常地图上使用的经纬度都为大地经纬度。

大地经度：参考椭球面上某点的大地子午面与本初子午面间的两面角。向东为正，向西为负。

大地纬度 ：参考椭球面上某点的法线与赤道平面的夹角。向北为正，向南为负。

大地高： 指某点沿法线方向到参考椭球面的距离。

只需要参考椭球体参数以及大地基准面就可以确定地理坐标系。

下面是Arcgis中对北京1954坐标系的说明。——WKID:4214

主要就是以下几个参数：

Prime Meridian（起始经度）

Datum（大地基准面）: D_Beijing_1954

Spheroid（参考椭球体）: Krasovsky_1940 （克拉索夫斯基椭球体）

西安-80地理坐标系。——WKID:4610

WGS-84地理坐标系。——WKID:4326

3.投影坐标系

在地球椭球面和平面之间建立点与点之间函数关系的数学方法，称为地图投影。

地球椭球表面是一种不可能展开的曲面，要把这样一个曲面表现到平面上，就会发生裂隙或褶皱。在投影面上，可运用经纬线的“拉伸”或“压缩”（通过数学手段）来加以避免，以便形成一幅完整的地图。地图投影的变形通常有：长度变形、面积变形和角度变形。在实际应用中，根据使用地图的目的，限定某种变形。

北京-54投影坐标系。——WKID:2435

国家2000投影坐标系。——WKID:4547

西安-80投影坐标系。——WKID:2383

WGS-84投影坐标系.。——WKID:3395

按变形性质分类：

等角投影：角度变形为零（Mercator）

等积投影：面积变形为零（Albers）

任意投影：长度、角度和面积都存在变形

其中，各种变形相互联系相互影响：等积与等角互斥，等积投影角度变形大，等角投影面积变形大。

从投影面类型划分：

横圆柱投影：投影面为横圆柱

圆锥投影：投影面为圆锥

方位投影：投影面为平面

从投影面与地球位置关系划分为：

正轴投影：投影面中心轴与地轴相互重合

斜轴投影：投影面中心轴与地轴斜向相交

横轴投影：投影面中心轴与地轴相互垂直

相切投影：投影面与椭球体相切

相割投影：投影面与椭球体相割

参考文献


展开全文
• 地理坐标和屏幕坐标相互转换，代码详细，使用
• 本包实现了WGS1984转本地坐标算法，对实现地理坐标转投影坐标及本地坐标系具有很好的参考价值，改包通过C#实现简单易懂。
• 你必须知道的地理坐标系和投影坐标系
仅供参考，欢迎指正。
1、基本概念
地理坐标系：为球面坐标。 参考平面地是椭球面，坐标单位：经纬度；
投影坐标系：为平面坐标。参考平面地是水平面，坐标单位：米、千米等；
地理坐标转换到投影坐标的过程可理解为投影。（投影：将不规则的地球曲面转换为平面）
2、地理坐标系
2.1 地球的三级逼近
2.1.1大地水准面
地球的自然表面有高山也有洼地，是崎岖不平的，我们要使用数学法则来描述他，就必须找到一个相对规则的数学面。
大地水准面是地球表面的第一级逼近。假设当海水处于完全静止的平衡状态时，从海平面延伸到所有大陆下部，而与地球重力方向处处正交的一个连续、闭合的曲面，这就是大地水准面。

2.1.2地球椭球体
大地水准面可以近似成一个规则成椭球体，但并不是完全规则，其形状接近一个扁率极小的椭圆绕短轴旋转所形成的规则椭球体，这个椭球体称为地球椭球体。它是地球的第二级逼近。

下面列举了一些常见椭球体的参数。我国1952年以前采用海福特椭球体，从1953年起采用克拉索夫斯基椭球体。 1978年我国决定采用新椭球体GRS（1975），并以此建立了我国新的、独立的大地坐标系，对应ArcGIS里面的Xian_1980椭球体。从1980年开始采用新椭球体GRS（1980），这个椭球体参数与ArcGIS中的CGCS2000椭球体相同。

2.1.3大地基准面
确定了一个规则的椭球表面以后，我们会发现还有一个问题，参考椭球体是对地球的抽象，因此其并不能去地球表面完全重合，在设置参考椭球体的时候必然会出现有的地方贴近的好（参考椭球体与地球表面位置接近），有地地方贴近的不好的问题，因此这里还需要一个大地基准面来控制参考椭球和地球的相对位置。 这是地球表面的第三级逼近。有以下两类基准面：
地心基准面：由卫星数据得到，使用地球的质心作为原点，使用最广泛的是 WGS 1984。
区域基准面：特定区域内与地球表面吻合，大地原点是参考椭球与大地水准面相切的点，例如Beijing54、Xian80。我们通常称谓的Beijing54、Xian80坐标系实际上指的是我国的两个大地基准面。
我们通常说的参心大地坐标系和地心大地坐标系的区别就在于此。
参心大地坐标系：指经过定位与定向后，地球椭球的中心不与地球质心重合而是接近地球质心。区域性大地坐标系。是我国基本测图和常规大地测量的基础。如Beijing54、Xian80。
地心大地坐标系：指经过定位与定向后，地球椭球的中心与地球质心重合。如CGCS2000、WGS84。
2.2地理坐标
地理坐标，就是用经线（子午线）、纬线、经度、纬度表示地面点位的球面坐标。
一般地理坐标可分为三种，天文经纬度，大地经纬度，地心经纬度。通常地图上使用的经纬度都为大地经纬度，所以这里我介绍一下大地经纬度，其他两种要想了解的话可以百度一下，其实区别不大。
大地经纬度：

大地经度：参考椭球面上某点的大地子午面与本初子午面间的两面角。东正西负。
大地纬度  ：参考椭球面上某点的法线与赤道平面的夹角。北正南负。
大地高： 指某点沿法线方向到参考椭球面的距离。
看到这里，地理坐标系的思路基本明确的了吧！只需要参考椭球体参数以及大地基准面就可以确定地理坐标系。下面是Arcgis中对北京1954坐标系的说明。

主要就是以下几个参数：
Prime Meridian（起始经度）
Datum（大地基准面）: D_Beijing_1954
Spheroid（参考椭球体）: Krasovsky_1940 （克拉索夫斯基椭球体）

3、投影坐标系
我们在选择坐标系的时候经常会发现以下情况：

这一大堆1954坐标系究竟是什么鬼,beijing1954不是地理坐标系吗？
为什么投影坐标系里也有？相信懵逼的不止我一个···
首先，投影坐标系的生成是以地理坐标系为基准的，所以每个投影坐标系前面都会挂有地理坐标系。而地理坐标系后面的一串乱七八糟的，则是投影参数！
比如  Beijing 1954 3 Degree GK Zone 39
意思是
3度分带法的北京54坐标系，中央经线在东117度的分带坐标，横坐标前加带号。
3.1投影
在地球椭球面和平面之间建立点与点之间函数关系的数学方法，称为地图投影。
地球椭球表面是一种不可能展开的曲面，要把这样一个曲面表现到平面上，就会发生裂隙或褶皱。在投影面上，可运用经纬线的“拉伸”或“压缩”（通过数学手段）来加以避免，以便形成一幅完整的地图。但不可避免会产生变形。
地图投影的变形通常有：长度变形、面积变形和角度变形。在实际应用中，根据使用地图的目的，限定某种变形。
根据不同的需要，我们会选择不同的投影组合！
按变形性质分类：
等角投影：角度变形为零（Mercator）
等积投影：面积变形为零（Albers）
任意投影：长度、角度和面积都存在变形
其中，各种变形相互联系相互影响：等积与等角互斥，等积投影角度变形大，等角投影面积变形大。
从投影面类型划分：

横圆柱投影：投影面为横圆柱
圆锥投影：投影面为圆锥
方位投影：投影面为平面
从投影面与地球位置关系划分为：

正轴投影：投影面中心轴与地轴相互重合
斜轴投影：投影面中心轴与地轴斜向相交
横轴投影：投影面中心轴与地轴相互垂直
相切投影：投影面与椭球体相切
相割投影：投影面与椭球体相割


投影参数：
标准线
概念：投影面与参考椭球的切线或割线。分为标准纬线与标准经线。
特点：没有变形，也称主比例尺。
中心线
概念：是指中央经线（原点经线）与中央纬线（原点纬线），用来定义图投影的中心或者原点。
特点：一般会有变形。

3.2我国常用投影
3.2.1高斯-克吕格投影
我国基本比例尺地形图（1：100万、1：50万、1：25万、1：10万、1：5万、1：2.5万、1：1万、1：5000）除1：100万以外均采用高斯-克吕格Gauss-Kruger投影（横轴等角切圆柱投影）为地理基础。
高斯克吕格投影的特点：
横轴等角切圆柱投影
– 离开中央子午线越远，变形越大
– 赤道是直线，离开赤道的纬线是弧线，凸向赤道
– 没有角度变形
– 长度和面积变形很小
北京54和西安80投影坐标系的投影方式
高斯投影特点:
– 中央子午线长度变形比为1
– 在同一条经线上，长度变形随纬度的降低而增大，在赤道处为最大
– 在同一条纬线上，长度变形随经差的增加而增大，且增大速度较快

我们经常会听到6°分带，3°分带的说法。其实并不是所有投影都有分带，从下面一张图就可以看出，分带是高斯克吕格投影自带的。

高斯－克吕格投影分带规定：该投影是国家基本比例尺地形图的数学基础，为控制变形，采用分带投影的方法，在比例尺1：2.5万—1：50万图上采用6°分带，对比例尺为1：1万及大于1：1万的图采用3°分带。
6°分带法：从格林威治零度经线起，每6°分为一个投影带，全球共分为60个投影带，东半球从东经0°—6°为第一带，中央经线为3°，依此类推，投影带号为1—30。其投影代号n和中央经线经度L0的计算公式为：L0=(6n—3)°；西半球投影带从180°回算到0°，编号为31—60，投影代号n和中央经线经度L0的计算公式为L0=360—(6n—3)°。
3°分带法：从东经1°30′起，每3°为一带，将全球划分为120个投影带，东经1°30′—4°30′，…178°30′—西经178°30′，…1°30′—东经1°30′。
东半球有60个投影带，编号1—60，各带中央经线计算公式：L0=3°n，中央经线为3°、6°…180°。西半球有60个投影带，编号1—60，各带中央经线计算公式：L0=360°—3°n，中央经线为西经177°、…3°、0°。
为了便于地形图的测量作业，在高斯-克吕格投影带内布置了平面直角坐标系统，具体方法是，规定中央经线为X轴，赤道为Y轴，中央经线与赤道交点为坐标原点，x值在北半球为正，南半球为负，y值在中央经线以东为正，中央经线以西为负。由于我国疆域均在北半球，x值均为正值，为了避免y值出现负值，规定各投影带的坐标纵轴均西移500km，中央经线上原横坐标值由0变为500km。为了方便带间点位的区分，可以在每个点位横坐标y值的百千米位数前加上所在带号。

3.2.2其他投影
1：100万地形图采用兰伯特Lambert投影（正轴等角割圆锥投影），其分幅原则与国际地理学会规定的全球统一使用的国际百万分之一地图投影保持一致。
海上小于50万的地形图多用墨卡托Mercator投影（正轴等角圆柱投影）。


展开全文
• 原始地理坐标是大地坐标 国内的GPS设备采集输出的坐标，应该是已经经过加偏的火星坐标吧 百度坐标，应该是在火星坐标的基础上二次加偏（或者是纠偏）的坐标 这样理解有问题嘛 把GPS传来的经纬度坐标经过百度的接口...
• dll中的方法及其参数有详细描述。 支持地理坐标向投影坐标的转换，以及投影坐标向地理坐标的转换。 主要输入参数为WKID。 如地理坐标系WGS1984的4326转成投影坐标系Beijing1954的2433
• 这是从网站上获取的地理坐标系和大地坐标系列表文件,文件格式为MHT
• dll中方法及其参数均有详细描述。 地理坐标，即经纬度，参数geoWkid=4326的WGS1984 投影坐标，即平面XY，参数prjWkid=2433的北京54
• 老师给我滴！！ARCGIS中坐标转换及地理坐标、投影坐标定义
• 包含打开地图，添加shp图层，保存地图的代码。关键是新建shp时用户自有选择地理坐标系和投影坐标系。坐标系的选择是通过prj文件来创建空间参考，因此本代码还带了arcgis坐标系的所有prj文件，方便选择坐标系。
• 支持地理坐标系变换，支持wgs84，北京54，西安80，csgc2000等坐标系下经纬度坐标和平面坐标互转。。
• 我使用GDAL库写了四个函数分别进行投影坐标与地理坐标（经纬度）之间的转换，投影坐标和图上坐标（行列号）之间的转换。有需要的朋友可以参考。 直接上代码吧，因为代码很简单（Python版本）。# -*- encoding: utf-...
我使用GDAL库写了四个函数分别进行投影坐标与地理坐标（经纬度）之间的转换，投影坐标和图上坐标（行列号）之间的转换。有需要的朋友可以参考。
直接上代码吧，因为代码很简单（Python版本）。

# -*- encoding: utf-8 -*-

from osgeo import gdal
from osgeo import osr
import numpy as np

def getSRSPair(dataset):
'''
获得给定数据的投影参考系和地理参考系
:param dataset: GDAL地理数据
:return: 投影参考系和地理参考系
'''
prosrs = osr.SpatialReference()
prosrs.ImportFromWkt(dataset.GetProjection())
geosrs = prosrs.CloneGeogCS()
return prosrs, geosrs

def geo2lonlat(dataset, x, y):
'''
将投影坐标转为经纬度坐标（具体的投影坐标系由给定数据确定）
:param dataset: GDAL地理数据
:param x: 投影坐标x
:param y: 投影坐标y
:return: 投影坐标(x, y)对应的经纬度坐标(lon, lat)
'''
prosrs, geosrs = getSRSPair(dataset)
ct = osr.CoordinateTransformation(prosrs, geosrs)
coords = ct.TransformPoint(x, y)
return coords[:2]

def lonlat2geo(dataset, lon, lat):
'''
将经纬度坐标转为投影坐标（具体的投影坐标系由给定数据确定）
:param dataset: GDAL地理数据
:param lon: 地理坐标lon经度
:param lat: 地理坐标lat纬度
:return: 经纬度坐标(lon, lat)对应的投影坐标
'''
prosrs, geosrs = getSRSPair(dataset)
ct = osr.CoordinateTransformation(geosrs, prosrs)
coords = ct.TransformPoint(lon, lat)
return coords[:2]

def imagexy2geo(dataset, row, col):
'''
根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换）
:param dataset: GDAL地理数据
:param row: 像素的行号
:param col: 像素的列号
:return: 行列号(row, col)对应的投影坐标或地理坐标(x, y)
'''
trans = dataset.GetGeoTransform()
px = trans[0] + col * trans[1] + row * trans[2]
py = trans[3] + col * trans[4] + row * trans[5]
return px, py

def geo2imagexy(dataset, x, y):
'''
根据GDAL的六 参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）
:param dataset: GDAL地理数据
:param x: 投影或地理坐标x
:param y: 投影或地理坐标y
:return: 影坐标或地理坐标(x, y)对应的影像图上行列号(row, col)
'''
trans = dataset.GetGeoTransform()
a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
b = np.array([x - trans[0], y - trans[3]])
return np.linalg.solve(a, b)  # 使用numpy的linalg.solve进行二元一次方程的求解

if __name__ == '__main__':
gdal.AllRegister()
dataset = gdal.Open(r"F:\2016\Data\Great Khingan\DEM\Projection\strm_6102_UTM.tif")
print('数据投影：')
print(dataset.GetProjection())
print('数据的大小（行，列）：')
print('(%s %s)' % (dataset.RasterYSize, dataset.RasterXSize))

x = 464201
y = 5818760
lon = 122.47242
lat = 52.51778
row = 2399
col = 3751

print('投影坐标 -> 经纬度：')
coords = geo2lonlat(dataset, x, y)
print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1]))
print('经纬度 -> 投影坐标：')
coords = lonlat2geo(dataset, lon, lat)
print('(%s, %s)->(%s, %s)' % (lon, lat, coords[0], coords[1]))

print('图上坐标 -> 投影坐标：')
coords = imagexy2geo(dataset, row, col)
print('(%s, %s)->(%s, %s)' % (row, col, coords[0], coords[1]))
print('投影坐标 -> 图上坐标：')
coords = geo2imagexy(dataset, x, y)
print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1]))


输出结果：

数据投影：
PROJCS["WGS_1984_UTM_Zone_51N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32651"]]
数据的大小（行，列）：
(7503 4799)
投影坐标 -> 经纬度：
(464201, 5818760)->(122.472422555, 52.5177753994)
经纬度 -> 投影坐标：
(122.47242, 52.51778)->(464200.830381, 5818760.513)
图上坐标 -> 投影坐标：
(2399, 3751)->(464163.754715, 5818797.73095)
投影坐标 -> 图上坐标：
(464201, 5818760)->(2399.49875769, 3751.50526134)

注：关于投影坐标和图上坐标转换的六参数模型可以参考我的另外一篇博文：经纬度坐标和投影坐标的转换，其实质就是一个仿射变换。

我们可以使用GDAL库自带的命令行工具（gdallocationinfo）进行检测：

其中参数-geoloc表示的后面给定坐标是投影坐标，-wgs84表示是WGS84参考系下的地理坐标（经纬度）。其输出是对应的图上坐标（行列号）。
具体参数可以使用gdallocationinfo –help查看。
展开全文
• 说明：椭球体、基准面构成了地理坐标系，即大地坐标系（经纬度）； 椭球体、基准面、投影构成了投影坐标系，即平面坐标系（米）。
• 深圳地铁站的地理坐标，包括地铁站名称，详细地址，经纬度信息
• 地心坐标系转换为地理坐标系，MATLAB程序，简单易懂，适合初学者。
• 全国每个城市的具体地理坐标，很精确。文件格式是.xml的。例如辽宁省北票市">41.800684#120.77073 # 是分割符。
• 该mat文件用于matlab中将经纬度坐标转换成地理坐标，按照北京54，第20度带进行坐标到平面的投影，调用projfwd函数即可。
• arcgis js地理坐标与屏幕坐标互转换
• 一般情况下，图层的坐标点由经纬度表示，单位为度，这是地理坐标系（地理坐标系是地球椭球体上的坐标，用经纬度表示）。但是当需要计算距离、面积等属性的时候，坐标点的单位必须是长度单位，这是投影坐标系（投影...
• echarts-all.js中的所有省市的地理坐标数组，以javascript数组的形式给出，可以直接引用
• //经纬度转墨卡托 -(CGPoint )lonLat2Mercator:(CGPoint ) lonLat { CGPoint mercator; double x = lonLat.x *20037508.34/180; double y = log(tan((90+lonLat.y)*M_PI/360))/(M_PI/180);...
• 用于经纬度地理坐标度、秒、分等各种形式相互之间的转换。
• 地理坐标系 3.1大地水准面 3.2 地球椭球体 3.3 大地基准面 四.投影坐标系 4.1 投影 4.2 投影带的计算 4.2.16°分带法 4.2.2 3°分带法 一.引言 地理信息系统（Geographic Information System或 Geo－...

...