基于MATLAB的等差分緯線多圓錐投影變換

摘 要：在MATLAB平臺下，綜合利用等差分緯線多圓錐投影的正反解變換實現了投影后矢量和柵格數據的輸出。輸出后的成圖不僅保證了數據的準確性，可靠性，并且可以供用戶直接在地理信息系統軟件ArcMap中讀取和編輯，方便用戶直接使用。

關鍵詞：等差分緯線多圓錐投影 MATLAB 投影變換 TFW

中圖分類號：P226 文獻標識碼：A 文章編號：1672-3791（2013）06（b）-0039-03

等差分緯線多圓錐投影是一種屬于任意性質的多圓錐投影，由我國地圖工作者于1963年根據我國同周圍的國家和地區的位置和形狀關系自主設計的一種投影。使用這種投影制作世界地圖時，往往選取東經150°作為中央經線，這樣不僅能完整顯示太平洋和各洲大陸，還能使我國版圖處于圖上的相對居中位置[1，2]。因此該投影在中國大陸出版的各種世界地圖專題圖中被廣泛應用，取得了較好的效果。

等差分緯線多圓錐投影的性質為：投影后的緯線（包括極點）為同軸圓弧，每個圓弧的圓心都位于中央經線上；中央經線的投影為一直線，其他經線的投影為曲線，并與中央經線對稱，且離中央經線越遠，其經線的間隔越成比例地遞減；極點投影長度為赤道投影長度的一半。

由于等差分緯線多圓錐投影沒有給出直接的投影公式，所以用戶在使用該投影時往往需要自己來設置一些參考點才能夠進行曲線擬合[3]，求解出該投影的正反解變換公式。

MATLAB是一種集數學計算功能、圖形化顯示功能等多種功能于一體的高級科學計算軟件，能夠高效快速地解決各種數學相關的科學問題，因此MATLAB被廣泛應用于不同領域[4]。本文將在MATLAB平臺下進行等差分緯線多圓錐投影的正反解變換，實現該投影下世界地圖的矢量和柵格數據輸出，以使這些數據可以直接在地理信息系統軟件ArcMap中讀取和編輯，方便用戶直接使用。

1 等差分緯線多圓錐投影變換在MATLAB中的實現

1.1 投影變換使用的地圖源數據

（1）世界地圖的.shp格式文件（見圖1）。

（2）像元行列比為1∶2的柵格世界地圖（8192×4096像素）（見圖2）。

1.2 投影變換使用的函數文件

以下均為在MATLAB軟件平臺下進行投影變換所用到的公式。這些公式可以在MATLAB軟件中各自存儲成一個.m類型的函數文件，以方便其他文件調用。

（1）正解變換公B635cr2INJ0j3S4ZldtYyrv2oU3pcYuQGWpJzFm57HQ=式：（xi，yi）=map_forwar d（long，lati）

式中，long為經度，范圍在-180～180之間；lati為緯度，范圍在-90-90之間；在經度值為long和緯度值為lati的條件下，xi，yi分別為經等差分緯線多圓錐投影正解變換[1，3，5，6]求得的x值和y值。若long和lati的值超過一定范圍時，xi和yi均返回NaN值。

（2）反解變換公式：（long，lati）=map_inverse（xi，yi）

式中各個變量的定義和前面第2節（1）相同，在x、y坐標值分別為xi和yi的條件下long和lati分別為經等差分緯線多圓錐投影反解變換[1，6]求得的經度值和緯度值。若xi和yi的值超過投影區域范圍時，long和lati均返回NaN值。

（3）經緯度轉像元行列號的公式：（x，y）=image_position（long，lati，image）

式中long和lati的定義和前面第2節（1）相同，image為讀入的圖像矩陣，在該函數的內部讀取該圖像的總行數imx和總列數imy。像元行列號x和y分別與long和lati的換算關系如下：

（1）

（2）

式中，[x]代表不超過x的最大正整數。

該式通過指定的經緯度轉換成柵格數據所對應的像元的行列號，把該行列號所在的像元值賦值給指定區域。因此該式主要用于柵格數據等差分緯線多圓錐投影的變換。

1.3 投影變換在MATLAB中的實現

（1）矢量數據的投影變換。

以前面所述的.shp矢量世界地圖為例，設使用的矢量世界地圖文件名為“world _map.shp”，讀取、轉換、顯示和輸出中央經線為東經150°的等差分緯線多圓錐投影世界地圖代碼如下[7]：

S=shaperead（'world_map.shp'）； %讀取shapefile文件。

Slength=length（S）； %讀取屬性長度。

for cou=1：Slength

xlength=length（S（cou，1）.X）；

for count=1：xlength

long=S（cou，1）.X（count）；lati=S（cou，1）.Y（count）； %讀取某一屬性的經緯度坐標。

if long<=-30 %此條件語句為輸出中央經線為東經150的投影時

long=long+210； %的經度變換。

else

long=long-150；

end

[xi，yi]=map_forward（long，lati）； %遍歷該Shapefile文件上的每一點，通過逐

S（cou，1）.X（count）=yi； %點變換，最終生成該投影下的地圖數據。

S（cou，1）.Y（count）=xi；

end

end

mapshow（S）； %顯示投影變換后的地圖數據。

shapewrite（S，'world_map_conversion.shp'）；%導出地圖數據。

在以上代碼中，第二行執行后的S變量實際上是一個數組，且數組中的每一個元素都是結構體，上面顯示了幾何類型、屬性名稱、X坐標和Y坐標等信息。通過嵌套循環語句逐個訪問了每一個點的經緯度坐標，并將變換后的X、Y坐標逐個賦值給S中原來的經緯度坐標。變換后還輸出并保存為Shapefile格式文件。變換后的結果如圖3所示。

（2）柵格數據的投影變換。

以前面所述的柵格世界地圖數據為例，設使用的柵格世界地圖文件名為“world_map.jpg”，輸出分辨率為5000×4000的中央經線為東經150°的等差分緯線多圓錐投影世界地圖。實現這一操作的代碼如下：

Mapc=uint8（ones（4000，5000，3））； %建立三維矩陣，圖像矩陣類型為uint8。

imagec=imread（'world_map.jpg'）； %讀取柵格世界地圖文件。

Mapc=Mapc*255； %使Mapc輸出的圖片背景為白色。

[p，q，r]=size（Mapc）； %讀取Mapc矩陣的尺寸。

for i=1：p

xi=（i-2000）*7000； %第6行和第8行：通過一系列變換，使xi和yi

for j=1：q %的范圍能包含投影范圍。

yi=（j-2500）*7000；

[long，lati]=map_inverse（xi，yi）； %投影反解變換。

if isnan（long）&&isnan（lati） %若long和lati的值為NaN，不進行顏色賦值，

continue； %并且直接跳入下一輪循環。

end

if long<=30&&long>=-180 %此條件語句為輸出中央經線為東經150的投影。

long=long+150；

else

if long<=180

long=long-210；

end

end

[row，column]=image_position（long，lati，imagec）； %讀取經緯度對應柵格世界

for k=1：r %地圖的fbd1f2ec6d0e10dd66a89dd291e25f47行列號。

Mapc（i，j，k）=imagec（row，column，k）； %把世界地圖的任意一點的

end %像元值賦值給投影地圖的

end %對應位置上。

end

image（Mapc）； %顯示地圖圖像。

該代碼首先是建立5000×4000的空白圖像，然后通過嵌套循環語句逐點掃描該空白圖像，把每一個像元的行列號（i，j）經一系列的變換轉化成該投影下的x、y坐標，利用反解變換求解對應的經緯度坐標，之后讀取經緯度坐標所對應的柵格世界地圖圖像所對應的行列號（row，column），把世界地圖中行列號為（row，column）的像元值賦值給5000×4000圖像的行列號為（i，j）的像元。當循環語句結束后，顯示投影變換后的圖像并保存。由于保存后的圖像已縮小7000倍，這樣如果在ArcMap中直接打開該圖像，會導致數據顯示不正確。但“GeoTiff”解決了這一瓶頸。tfw文件是TIFF文件坐標信息的文本文件，ArcInfo、MicroStation、AutoCAD等軟件均支持該格式的坐標信息文件。此文件定義了柵格圖像素坐標與實際大地坐標的仿射關系[8，9]。一個.tfw格式的文件構成如下[10]：

第1行：A：地理x坐標中的像元分辨率

第2行：D：y軸旋轉系數

第3行：B：x軸旋轉系數

第4行：E：地理y坐標中的像元分辨率

第5行：C：左上角第1個像元中心的x坐標值

第6行：F：左上角第1個像元中心的y坐標值

各個參數的圖解如圖4所示。

由此保存分辨率為5000×4000的投影后的世界地圖為TIF格式，用記事本編寫對應的tfw格式文件。最后在ArcMap中加載這個世界地圖，如圖5所示。

2 結語

本文依托MATLAB軟件在圖形處理方面的優勢，綜合利用等差分緯線多圓錐投影的正反解變換公式，還利用了像素長寬比為2∶1的世界地圖原始數據中經緯度與行列號的關系，實現了矢量數據和柵格數據的等差分緯線多圓錐投影變換。此外本文以桌面GIS軟件ArcMap為例，使用了.tfw格式文件讓輸出后柵格數據等差分緯線多圓錐投影世界地圖準確顯示。本文中的投影變換呈現出如下特點。

（1）輸出后的矢量和柵格地圖數據可以供用戶在ArcMap、MapInfo等GIS軟件中被直接加載，方便用戶進行編輯。

（2）輸出后的矢量和柵格地圖數據具有準確性、可靠性。與文獻[3]相比，整個變換過程無需人工干預，節省了人工計算配準點的繁重勞動，減少了柵格地圖數據地理配準中可能出現的誤差。

參考文獻

[1]董曼，李勝樂.世界地圖等差分緯線多圓錐投影的正反解變換[J].大地測量與地球動力學，2008，28（2）：95-99.

[2]祝國瑞.地圖學[M].武漢：武漢大學出版社，2004.

[3]葉遠智.柵格數據的等差分緯線多圓錐投影轉換[J].測繪通報，2012（9）：68-70.

[4]劉正君.MATLAB科學計算寶典[M].北京：電子工業出版社，2012.

[5]胡毓鉅.地圖投影[M].北京：測繪出版社，1981.

[6]楊啟和.地圖投影變換原理和方法[M].北京：解放軍出版社，1989.

[7]蘇金明，王永利.MATLAB圖形圖像[M].北京：電子工業出版社，2005.

[8]余建軍，蔡以雷，洪景峰.一種利用數字線劃圖快速生產數字柵格地圖的方法[J].測繪通報，2009（7）：51-54.

[9]王憲民.ArcGIS在柵格數據處理中的應用初探[J].測繪技術裝備，2005，1（7）： 32-34.

[10]吳銘杰.基于AutoCAD的TIFF world（TFW）圖像的應用研究[A].史照良.第十四屆華東六省一市測繪學會學術交流會論文集[C].南京：《現代測繪》編輯部，2012：157-158.