gdal坐標轉換總結(轉換）

轉自https://blog.csdn.net/qq_32657025/article/details/80176520

首先，在進行坐標轉換之前，有必要先了解一下有關坐標系的幾個基本概念。

地理坐標系（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 =@null +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 +nadgrids=@null +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中的顯示如下圖：

理論上講，”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工具為影響重新定義投影信息：

• gdal_translate:會將重定義投影后的影像存在新的影像文件中；
• gdal_edit.py：可以直接修改指定影像文件的投影信息。

參考文獻和測試數據

Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
Revisiting Web Mercator
Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
embeding esri compatible arcgis projection using gdal

編譯好的各版本gdal庫下載

gdal官方幫助文檔
EPSG查詢官網

《基於Proj_4的空間坐標轉換》蓋森

測試數據和參考文獻下載鏈接：https://pan.baidu.com/s/1RE3e-9Qce2QRDHYcoEvbew 密碼：yg1d