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

關(guān)鍵詞：等差分緯線多圓錐投影 MATLAB 投影變換 TFW

中圖分類號(hào)：P226 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1672-3791（2013）06（b）-0039-03

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

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

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

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

1 等差分緯線多圓錐投影變換在MATLAB中的實(shí)現(xiàn)

1.1 投影變換使用的地圖源數(shù)據(jù)

（1）世界地圖的.shp格式文件（見(jiàn)圖1）。

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

1.2 投影變換使用的函數(shù)文件

以下均為在MATLAB軟件平臺(tái)下進(jìn)行投影變換所用到的公式。這些公式可以在MATLAB軟件中各自存儲(chǔ)成一個(gè).m類型的函數(shù)文件，以方便其他文件調(diào)用。

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

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

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

式中各個(gè)變量的定義和前面第2節(jié)（1）相同，在x、y坐標(biāo)值分別為xi和yi的條件下long和lati分別為經(jīng)等差分緯線多圓錐投影反解變換[1，6]求得的經(jīng)度值和緯度值。若xi和yi的值超過(guò)投影區(qū)域范圍時(shí)，long和lati均返回NaN值。

（3）經(jīng)緯度轉(zhuǎn)像元行列號(hào)的公式：（x，y）=image_position（long，lati，image）

式中l(wèi)ong和lati的定義和前面第2節(jié)（1）相同，image為讀入的圖像矩陣，在該函數(shù)的內(nèi)部讀取該圖像的總行數(shù)imx和總列數(shù)imy。像元行列號(hào)x和y分別與long和lati的換算關(guān)系如下：

（1）

（2）

式中，[x]代表不超過(guò)x的最大正整數(shù)。

該式通過(guò)指定的經(jīng)緯度轉(zhuǎn)換成柵格數(shù)據(jù)所對(duì)應(yīng)的像元的行列號(hào)，把該行列號(hào)所在的像元值賦值給指定區(qū)域。因此該式主要用于柵格數(shù)據(jù)等差分緯線多圓錐投影的變換。

1.3 投影變換在MATLAB中的實(shí)現(xiàn)

（1）矢量數(shù)據(jù)的投影變換。

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

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

Slength=length（S）； %讀取屬性長(zhǎng)度。

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）； %讀取某一屬性的經(jīng)緯度坐標(biāo)。

if long<=-30 %此條件語(yǔ)句為輸出中央經(jīng)線為東經(jīng)150的投影時(shí)

long=long+210； %的經(jīng)度變換。

else

long=long-150；

end

[xi，yi]=map_forward（long，lati）； %遍歷該Shapefile文件上的每一點(diǎn)，通過(guò)逐

S（cou，1）.X（count）=yi； %點(diǎn)變換，最終生成該投影下的地圖數(shù)據(jù)。

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

end

end

mapshow（S）； %顯示投影變換后的地圖數(shù)據(jù)。

shapewrite（S，'world_map_conversion.shp'）；%導(dǎo)出地圖數(shù)據(jù)。

在以上代碼中，第二行執(zhí)行后的S變量實(shí)際上是一個(gè)數(shù)組，且數(shù)組中的每一個(gè)元素都是結(jié)構(gòu)體，上面顯示了幾何類型、屬性名稱、X坐標(biāo)和Y坐標(biāo)等信息。通過(guò)嵌套循環(huán)語(yǔ)句逐個(gè)訪問(wèn)了每一個(gè)點(diǎn)的經(jīng)緯度坐標(biāo)，并將變換后的X、Y坐標(biāo)逐個(gè)賦值給S中原來(lái)的經(jīng)緯度坐標(biāo)。變換后還輸出并保存為Shapefile格式文件。變換后的結(jié)果如圖3所示。

（2）柵格數(shù)據(jù)的投影變換。

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

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行：通過(guò)一系列變換，使xi和yi

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

yi=（j-2500）*7000；

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

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

continue； %并且直接跳入下一輪循環(huán)。

end

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

long=long+150；

else

if long<=180

long=long-210；

end

end

[row，column]=image_position（long，lati，imagec）； %讀取經(jīng)緯度對(duì)應(yīng)柵格世界

for k=1：r %地圖的fbd1f2ec6d0e10dd66a89dd291e25f47行列號(hào)。

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

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

end %對(duì)應(yīng)位置上。

end

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

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

第1行：A：地理x坐標(biāo)中的像元分辨率

第2行：D：y軸旋轉(zhuǎn)系數(shù)

第3行：B：x軸旋轉(zhuǎn)系數(shù)

第4行：E：地理y坐標(biāo)中的像元分辨率

第5行：C：左上角第1個(gè)像元中心的x坐標(biāo)值

第6行：F：左上角第1個(gè)像元中心的y坐標(biāo)值

各個(gè)參數(shù)的圖解如圖4所示。

由此保存分辨率為5000×4000的投影后的世界地圖為T(mén)IF格式，用記事本編寫(xiě)對(duì)應(yīng)的tfw格式文件。最后在ArcMap中加載這個(gè)世界地圖，如圖5所示。

2 結(jié)語(yǔ)

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

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

（2）輸出后的矢量和柵格地圖數(shù)據(jù)具有準(zhǔn)確性、可靠性。與文獻(xiàn)[3]相比，整個(gè)變換過(guò)程無(wú)需人工干預(yù)，節(jié)省了人工計(jì)算配準(zhǔn)點(diǎn)的繁重勞動(dòng)，減少了柵格地圖數(shù)據(jù)地理配準(zhǔn)中可能出現(xiàn)的誤差。

參考文獻(xiàn)

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

[2]祝國(guó)瑞.地圖學(xué)[M].武漢：武漢大學(xué)出版社，2004.

[3]葉遠(yuǎn)智.柵格數(shù)據(jù)的等差分緯線多圓錐投影轉(zhuǎn)換[J].測(cè)繪通報(bào)，2012（9）：68-70.

[4]劉正君.MATLAB科學(xué)計(jì)算寶典[M].北京：電子工業(yè)出版社，2012.

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

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

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

[8]余建軍，蔡以雷，洪景峰.一種利用數(shù)字線劃圖快速生產(chǎn)數(shù)字柵格地圖的方法[J].測(cè)繪通報(bào)，2009（7）：51-54.

[9]王憲民.ArcGIS在柵格數(shù)據(jù)處理中的應(yīng)用初探[J].測(cè)繪技術(shù)裝備，2005，1（7）： 32-34.

[10]吳銘杰.基于AutoCAD的TIFF world（TFW）圖像的應(yīng)用研究[A].史照良.第十四屆華東六省一市測(cè)繪學(xué)會(huì)學(xué)術(shù)交流會(huì)論文集[C].南京：《現(xiàn)代測(cè)繪》編輯部，2012：157-158.