1. 程式人生 > >gdal座標轉換總結

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 +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