gdal坐标转换总结(转换）

首先，在进行坐标转换之前，有必要先了解一下有关坐标系的几个基本概念。

地理坐标系（Geographic Coordinate Systems）

地理坐标系是一个球面的坐标系统，以经纬度为单位，它由椭球体和大地基准面两个部分组成。

椭球体（spheroid）

我们要将地理信息以球面坐标系的方式表达，首先需要找到一个可以量化计算的椭球体。一个椭球体的确定需要以下参数：长半轴、短半轴、偏心率，其中偏心率可根据长短半轴计算得到。

例如，WGS84椭球的参数如下：

Spheroid(椭球名):"WGS_84";
Semimajor Axis(长半轴):6378137
Semimajor Axis(长半轴):6356752.3142
Inverse Flattening(扁率):1/298.2572236

大地基准面(datum)

有了椭球体以后，还需要一个大地基准面将这个椭球定位。

大地基准面（Geodetic datum），设计为最密合部份或全部大地水准面的数学模式。它由椭球体本身及椭球体和地表上一点(原点)之间的关系来定义。此关系能以 6个量来定义，通常是大地纬度、大地经度、原点高度、原点垂线偏差之两分量及原点至某点的大地方位角。

同一个椭球面，不同的地区由于关心的位置不同，需要最大限度的贴合自己的那一部分，因而大地基准面就会不同。

有了Spheroid和Datum两个基本条件，便可以确定一个地理坐标系统。

投影坐标系

将球面坐标转化为平面坐标的过程称为投影。因此，投影坐标系实质上是在地理坐标系的基础上通过投影得到的。投影坐标系其单位通常为m。

例如我国常用的高斯-克吕格投影,其通常是按6度和3度分带投影，1:2.5万－1:50万比例尺地形图采用经差6度分带，1:1万比例尺的地形图采用经差3度分带。具体分带法是：6度分带从本初子午线开始，按经差6度为一个投影带自西向东划分，全球共分60个投影带，带号分别为1－60；3度投影带是从东经1度30秒经线开始，按经差3度为一个投影带自西向东划分，全球共分120个投影带。为了便于地形图的测量作业，在高斯-克吕格投影带内布置了平面直角坐标系统，具体方法是，规定*经线为X轴，赤道为Y轴，*经线与赤道交点为坐标原点，x值在北半球为正，南半球为负，y值在*经线以东为正，*经线以西为负。由于我国疆域均在北半球，x值均为正值，为了避免y值出现负值，规定各投影带的坐标纵轴均西移500km，*经线上原横坐标值由0变为500km。为了方便带间点位的区分，可以在每个点位横坐标y值的百千米位数前加上所在带号，如20带内A点的坐标可以表示为YA=20 745 921.8m。

ArcGIS中常用的投影坐标系命名方式

Beijing 1954 3 Degree GK CM 117E

北京54坐标系 3度分带投影高斯克吕格 *经线东经117度横坐标前加带号
Beijing 1954 3 Degree GK Zone 25

北京54坐标系 3度分带投影高斯克吕格带号25 横坐标前加带号
Beijing 1954 GK Zone 13

北京54坐标系 6度分带投影高斯克吕格带号13 横坐标前加带号
Beijing 1954 GK Zone 13N

北京54坐标系 6度分带投影高斯克吕格带号13 横坐标前不加带号

ArcGIS动态投影

ArcMap中的Data的空间参考默认为第一个加载到当前工作区的那个文件的坐标系统，后加入的数据，如果和当前工作区坐标系统不相同，则ArcMap会自动做投影变换，把后加入的数据投影变换到当前坐标系统下显示。但此时数据文件所存储的数据并没有改变，只是显示形态上的变化。因此叫动态投影。表现这一点最明显的例子就是，在Export Data时，会让你选择是按this layer’s source data（数据源的坐标系统导出），还是按照the Data （当前数据框架的坐标系统）导出数据。

坐标系的通用描述方法

WKT

WKT，全程 well know text，是一种文本标记语言，用于表示矢量几何对象、空间参照系统及空间参照系统之间的转换，该格式由开放地理空间联盟(OGC)制定。

例如，WGS84地理坐标系的OGC WTK定义如下：

GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

ESRI使用的WTK描述与OGC标准的WTK存在区别，因此，有时候会需要用到ESRI WTK，例如，WGS84地理坐标系的ESRI WTK如下：

GEOGCS["GCS_WGS_1984",
    DATUM["D_WGS_1984",
        SPHEROID["WGS_1984",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]

在arcgis中，.prj文件中存储的即是要素的坐标系的WKT文本表示，可用记事本打开查阅。

EPSG

EPSG的英文全称是European Petroleum Survey Group，中文名称为欧洲石油调查组织，这个组织成立于1986年，2005年并入IOGP(International Association of Oil & Gas Producers)，中文名称为国际油气生产者协会。

EPSG大地测量参数数据集是由EPSG组织定义的坐标参考系统和坐标转换的结构化数据集，现由IOGP地理信息委员会的大地测量小组委员会维护。不同的EPSG编号能代表不同的坐标系统、椭球体、大地水准面等等。例如：

EPSG编号	类型	代表
EPSG:4326	Geodetic coordinate system(地理坐标系)	WGS 84
EPSG:2437	Projected coordinate system(投影坐标系)	Beijing 1954 / 3-degree Gauss-Kruger CM 120E

可以通过网址:epsg.io查询不同坐标系的EPSG代号，或查询EPSG代号的含义，还可在该页面查询EPSG坐标系的WKT文本或者Proj.4表示。

Proj.4

Proj.4是开源GIS著名的地图投影库。Proj.4使用参数的方式来描述坐标系统，基本格式为：“+param=value”。

例如，WGS84地理坐标系的Proj.4表示方法为：

"+proj=longlat +datum=WGS84 +no_defs"

其中，“+proj”表示投影类型代字，“longlat”代表该坐标系为地理坐标系；“+datum”表示大地基准面代字。

GDAL

GDAL（geospatial data abstraction libarary），是一个开源的栅格空间数据转换库，它拥有一系列命令行工具用于数据转换和处理。OGR是GDAL库的一部分，提供对矢量数据的支持。

gdalinfo

用于列举栅格数据集的相关信息。可以获取到的信息包括：

栅格数据格式
栅格大小（行列值）
坐标系统（wkt/proj.4）
角坐标、中心坐标
文件元数据
波段信息
……

命令行命令格式示例：

//gdalinfo [栅格文件名]
gdalinfo L14.tif

输出示例：

Driver: GTiff/GeoTIFF
Files: L14.tif
Size is 22016, 15360
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids
[email protected] +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (12621282.110448303000000,3656747.433162831700000)
Pixel Size = (9.554628535647032,-9.554628535647032)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (12621282.110, 3656747.433) (113d22'44.06"E, 31d11' 4.59"N)
Lower Left  (12621282.110, 3509988.339) (113d22'44.06"E, 30d 3' 0.28"N)
Upper Right (12831636.812, 3656747.433) (115d16' 6.80"E, 31d11' 4.59"N)
Lower Right (12831636.812, 3509988.339) (115d16' 6.80"E, 30d 3' 0.28"N)
Center      (12726459.461, 3583367.886) (114d19'25.43"E, 30d37' 8.42"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue

gdal_translate

转换栅格数据的格式。在准换过程中可以执行数据设置、重采样等操作。例如，可以利用gdal_translate工具重新定义栅格数据的投影（与arcgis中的define project功能类似）。重新定义投影之后的栅格存储在新的文件中，源文件的投影信息不变。
如下命令，将L14.tif影响的投影信息重新定义为testPolygon.prj文件定义的投影，并将重新投影后的影像重新存储为L14_rePro.tif。

gdal_translate -a_srs testPolygon.prj L14.tif L14_rePro.tif

gdalsrsinfo

以指定的格式列举给定的空间参考系的相关信息。

支持的输入格式包括：

GDAL/ORG支持的，且包含空间参考系信息影响数据集的文件名；
GDAL/ORG常用的空间参考系规范，包括EPSG
、PROJ.4、WKT等。

支持的输出格式包括（-o out_type 的可选参数）：

default：默认为prj.4 和wkt格式
all：所有可用的格式
wkt_all: 所有可用的wkt格式
proj4: PROJ.4字符串格式
wkt：OGC定义的WKT模式（完整版）
wkt_simple: OGC定义的WKT模式（精简版）
wkt_noct： OGC定义的WKT模式 (不带 OGC CT params)
wkt_esri: esri定义的WKT格式
mapinfo: mapinfo风格的坐标系统格式
xml: xml格式

例如，可以分别查看EPSG:3857投影的的WKT表示方法和ESRI wkt表示方法，输入命令：

gdalsrsinfo -o wkt EPSG:3857
gdalsrsinfo -o wkt_esri EPSG:3857

wkt输出结果为：

PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m [email protected] +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

wkt_esri输出结果为：

PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

由上述结果可知，ESRI在表示web墨卡托投影时的WKT与OGC标准的WKT存在一定的差别，主要在于wkt_esri多了几个描述参数。

gdalwarp

gdalwarp是一个用于图像镶嵌，重投影和变形的工具，它可以将影像重投影到任何受支持的投影。
该工具的重要参数包括：

-s_srs srs_def：源空间参考系。任意可以被 OGRSpatialReference.SetFromUserInput()函数调用支持的坐标系，都可以作为参数。包括EPSG定义的投影坐标系（PCS)和地理坐标系(PCS)，PROJ.4定义的投影声明，或者包含投影信息的 .prj文件的文件名。
-t_srs srs_def：目标空间参考系。可用的参数与原空间参考系一致。
srcfile:源文件
dstfile:投影后的目标文件

例如，”EPSG:4326”代表WGS84地理坐标系，”EPSG:4214”代表beijing54地理坐标系，用gdalwrap实现二者之间的转换如下：


gdalwarp -s_srs "EPSG:4326" -t_srs "EPSG:4214" dem_wgs84.tif dem_wgs84_to_beijing54.tif

epsg:3857投影的geotiff影像在arcgis10.4以下版本中显示错位的问题

在使用gdal为栅格影像进行重新定位的过程中，我发现，当geotiff格式的影像的投影格式为EPSG:3857(OGC WKT)时，在Arcgis10.3中打开时会发生错位。

经查阅，GDAL目前在GeoTIFF中编码EPSG:3857的方式使得这种GeoTIFF在ArcGIS中打开时以非常显着的方式错位，因为ArcGIS会将其解释为带有WGS84的椭球定义的Mercator_1SP，而不是球形公式。因此，出现错位的原因是因为arcgis10.3及以下版本不能正确解析非ESRI格式的EPSG:3857(“WGS 84 / Pseudo-Mercator”）投影影像。而并非gdal进行投影转换时出现错误,而是原影像在arcgis中的显示错误。

EPSG:3857投影的OGC WKT格式为：

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1],
    AXIS["X",EAST],
    AXIS["Y",NORTH]]

EPSG:3857投影的ESRI WKT格式为：

PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",
  GEOGCS["GCS_WGS_1984",
    DATUM["D_WGS_1984",
      SPHEROID["WGS_1984",6378137.0,298.257223563]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]],
  PROJECTION["Mercator_Auxiliary_Sphere"],
  PARAMETER["False_Easting",0.0],
  PARAMETER["False_Northing",0.0],
  PARAMETER["Central_Meridian",0.0],
  PARAMETER["Standard_Parallel_1",0.0],
  PARAMETER["Auxiliary_Sphere_Type",0.0],
  UNIT["Meter",1.0]]

由上述描述可以看出，ESRI WKT格式的web墨卡托投影较之OGC格式的WKT多了几个参数，正式这几个参数，告诉arcgis如何正确的显示web墨卡托投影。

为了更加准确的验证这一点，我们再arcgis中新建了一个多边形图层，将并将该图层的SRS设置为arcgis内置的EPSG:3857投影，即”WGS_1984_Web_Mercator_Auxiliary_Sphere”投影，并使新建多边形的范围与L14元数据中提供的栅格的范围一致。两者在arcgis中的显示如下图：
gdal坐标转换总结(转换）

理论上讲，”WGS 84 / Pseudo-Mercator”投影和”WGS_1984_Web_Mercator_Auxiliary_Sphere”投影的EPSG代号均为3857，arcgis应该按照同一方式显示。但事实上，arcgis10.3对非esri格式的3847投影，并没有按照伪墨卡托投影正确显示。

解决方法：
arcgis10.4以后的版本针对这一问题进行了解决，在arcgis10.4中加载”WGS 84 / Pseudo-Mercator”投影的影像时，会自动转换为ESRI格式的伪墨卡托投影，即”WGS_1984_Web_Mercator_Auxiliary_Sphere”投影。

如果希望在arcgis10.3及以下版本中正确显示影像，需要修改影像的投影信息，修改为Arcgis能支持的伪墨卡托投影格式。

第一种修改方式为，在Arcgis里面对L14.tif的投影信息进行重新定义，注意，arcgis中 project工具和define project工具的区别在于

project

Project工具对源文件的x-y坐标起作用，可将其转换至不同的坐标系统，生成新的要素类或栅格影像，同时不改变原有要素类或栅格。新文件不仅具有新的坐标系统，而且还具有不同的坐标系统标注。若需将有坐标系统的图层转换为不同坐标系统，可以使用Project工具。
define Project

Define Projection工具只改变要素类或者栅格的坐标系统标注，而不会影响其内部坐标，只适合用于具有未知坐标系统的数据集，或者因标注错误而需要更正的数据集。

我们遇到的问题是栅格文件的坐标系统标记方式不被arcgis所支持，因此，只需要将该文件的坐标系重新定义为arcgis支持的同一坐标系即可。因此，使用arcgis将L14.tif的坐标系重新定义为”WGS_1984_Web_Mercator_Auxiliary_Sphere”时，影像便能在arcgis中正确显示。如下图，对影像进行define project操作后，用于测试的多边形已经能够与影像重合。
gdal坐标转换总结(转换）