1. 程式人生 > >經緯度定義、經緯度格式、GDAL中地理座標轉換及地理座標螢幕顯示

經緯度定義、經緯度格式、GDAL中地理座標轉換及地理座標螢幕顯示

關於經緯度的相關定義:

經線:

        經線也稱子午線,是人類為度量方便而假設出來的輔助線,定義為地球表面連線南北兩極的大圓線上的半圓弧。任兩根經線的長度相等,相交於南北兩極點。每一根經線都有其相對應的數值,稱為經度。經線指示南北方向。

本初子午線 :

       本初子午線又稱“首子午線”或“零子午線”,是計算東西經度的起點,以經過英國倫敦東南格林尼治的經線為本初子午線,作為計算地理的起點和世界標準“時區”的起點。本初子午線以東為東經,以西為西經,全球經度測量均以本初子午線與赤道的交點E點作為經度原點。

 子午線命名的由來:

      某一天體視運動軌跡中,同一子午線上的各點該天體在上中天(午)與下中天(子)出現的時刻相同。

緯線:

    緯線是人類為度量方便而假設出來的輔助線,定義為地球表面某點隨地球自轉所形成的軌跡。緯線平行於赤道,經線垂直於赤道,任何一根緯線都是圓形而且兩兩平行。緯線的長度是赤道的周長乘以緯線的緯度的餘弦,所以赤道最長,離赤道越遠的緯線,周長越短,到了兩極就縮為0。緯線指示東西方向。

北迴歸線(23°26' 22" N)

      太陽在北半球能夠直射到的離赤道最遠的位置,其緯度值為黃赤交角,是一條緯線,大約在北緯23度26分的地方。

赤道(0°N)

     通過地球中心劃一個與地軸成直角相交的平面,在地球表面相應出現一個和地球的極距離相等的假想圓圈。赤道的緯度是0°。是地球表面的點隨地球自轉產生的軌跡中周長最長的圓周線,赤道半徑 6378.137Km ;兩極半徑 6359.752Km;平均半徑 6371.012Km ;赤道周長 40075.7Km。如果把地球看做一個絕對的球體的話,赤道距離南北兩極相等,是一個大圓。它把地球分為南北兩半球,其以北是北半球,以南是南半球,是劃分緯度的基線,是南北緯線的起點。

南迴歸線(23° 26' 22" S)

       太對應於北迴歸線,陽在地球上的直射點在一年內到達的最南點所在的緯線。

經度:

     是地球上一個地點離一根被稱為本初子午線的南北方向走線以東或以西的度數。它是通過某地的經線面與本初子午面所成的二面角,在本初子午線以東的經度叫東經,在本初子午線以西的叫西經。東經用“E”表示,西經用“W”表示。比如:E116.33533°,W116.33533°.經度的每一度被分為60角分,每一分被分為60秒。也可以用小數點位表示。有時西經也被寫成負數,如:-116.33533°。

緯度

      是指某點與地球球心的連線和地球赤道面所成的線面角,其數值在0至90度之間。位於赤道以北的點的緯度叫北緯,記為N;位於赤道以南的點的緯度稱南緯,記為S。

緯度數值在0至30度之間的地區稱為低緯度地區;緯度數值在30至60度之間的地區稱為中緯度地區;緯度數值在60至90度之間的地區稱為高緯度地區。赤道是0°緯線,北緯度的最大值為90°,即北極點;南緯度的最大值為90°,即南極點。

地面法線

     垂直於某地參考橢球體表面的線。

經緯度座標系統:

        是經度與緯度的合稱組成一個座標系統。又稱為地理座標系統,它是一種利用三度空間的球面來定義地球上的空間的球面座標系統,能夠標示地球上的任何一個位置。

一個經度和一個緯度一起確定地球上一個地點的精確位置。一個完整的圖面地理座標系統由一個完整的參考橢球體系和一個完整的投影體系組成。

天文座標系統和大地座標系統

       由經線和緯線構成兩組互相正交的曲線座標網叫地理座標網。由經緯度構成的地理座標系統又叫地理座標系。按照是否原點是否在地球上,座標系可分為天文座標和大地地理座標。天文座標是用天文測量方法確定的,大地地理座標是用大地測量方法確定的。我們在地球橢球面上所用的地理座標系屬於大地地理座標系,簡稱大地座標系。

經緯度網

       大地座標系是以以參考橢球面為基準面的座標,地面點P的位置用大地經度L、大地緯度B和大地高H表示。當點在參考橢球面上時,僅用大地經度和大地緯度表示。雖然如此,但是,大地座標系與經緯度網並不一致,大地座標系以球面的形式描述了真實的地球座標,而經緯度網則是在製圖時,為滿足平面作圖的需要(即地圖投影),參照一定的投影公式製作成平面地圖,並在圖上套上經緯線。常用作世界地圖的投影有墨卡託投影、高爾投影、摩爾威特投影、等差分緯線多圓錐投影、格靈頓投影、桑森投影、烏爾馬耶夫投影等。而橫軸墨卡託投影即為高斯-克呂格投影,一般用於我國的地理製圖。
 我國的高斯-克呂格投影體系

       它是一種橫軸、橢圓柱面、等角投影。其投影過程可簡述如下:橢圓柱面與地球橢球在某一子午圈上相切,這條子午圈叫做投影的中央子午線,又稱軸子午線,它也是高斯投影后的平面直角座標系的縱軸(一般定義為x軸);地球的赤道面與橢圓柱面相交成一條直線,這條直線與中央子午線正交,它是平面直角座標系的橫軸(y軸);把橢圓柱面展開,就得出以(x,y)為座標的平面直角座標系。按一定經差將地球橢球面劃分成若干投影帶,這是高斯投影中限制長度變形的最有效方法。分帶時既要控制長度變形使其不大於測圖誤差,又要使帶數不致過多以減少換帶計算工作,據此原則將地球橢球面沿子午線劃分成經差相等的瓜瓣形地帶,以便分帶投影。通常按經差6度或3度分為六度帶或三度帶。六度帶自0度子午線起每隔經差6度自西向東分帶,帶號依次編為第1、2…60帶。三度帶是在六度帶的基礎上分成的,它的中央子午線與六度帶的中央子午線和分帶子午線重合,即自1.5度子午線起每隔經差3度自西向東分帶,帶號依次編為三度帶第 1、2…120帶。我國的經度範圍西起73°東至135°,可分成六度帶十一個,各帶中央經線依次為75°、81°、87°、……、117°、123°、129°、135°,或三度帶二十二個。

經緯度格式

經緯度格式分為三種:度、度-分、度-份-秒

1.   ddd.ddddd °【度 . 度 格式】的十進位制小數部分(5位)

2.   ddd°mm.mmm’ 【度 . 分 . 分 格式】的十進位制小數部分(3位)

3.  ddd°mm’ss’’ 【度 . 分 . 秒 格式】

轉換公式

度分轉換:
將度分單位資料轉換為度單位資料
度=度+分/60
例如:
經度 = 116°20.12’
緯度 = 39°12.34’
經度 = 116 + 20.12 / 60 = 116.33533°
緯度 = 39 + 12.34 / 60 = 39.20567°

度分秒轉換:
將度分秒單位資料轉換為度單位資料
度 = 度 + 分 / 60 + 秒 / 60 / 60
例如:
經度 = 116°20’43”
緯度 = 39°12’37”
經度 = 116 + 20 / 60 + 43 / 60 / 60 = 116.34528°
緯度 = 39 + 12 / 60 + 37 / 60 / 60 = 39.21028°

在實際的處理中,計算相應像元點的地理座標:

GDAL資料集有兩種模式描述柵格位置(用點/線座標系)以及地理參考座標系之間的關係:首要的也是最普遍的是使用仿射轉換,另一種則是GCPs(多控制點定位方式)

仿射變換

1用GDAL獲取四個角點的地理座標

GDALDataset::::GetGeoTransform()

宣告一個儲存變數:

double        geoTransform[6];

geoTransform[0]是左上角像元的東座標;

geoTransform[3]是左上角像元的北座標;

geoTransform[1]是影像寬度上的解析度;

geoTransform[5]是影像高度上的解析度;

geoTransform[2]是旋轉, 0表示上面為北方;

geoTransform[4]是旋轉, 0表示上面為北方;

2相應的放射變換公式:

Xp = geoTransform [0] +Xpixel*geoTransform [1]+Yline*geoTransform [2];

Yp = geoTransform [3] + Xpixel*geoTransform [4] + YlineL*geoTransform [5];

其中Xpixel、YlineL分別地物表示單純地影像上影象遍歷時(x,y)座標表示,點/線座標系是從左上角(0,0)點到右下角,也就是座標軸從左到右增長,從上到下增長的座標系。

2 GCPs控制點座標變換:

int nGCPs=poDataset->GetGCPCount();   //獲得控制點數目

const GDAL_GCP* pGCPs = poDataset->GetGCPs(); //獲得GCP控制點,是一個欄位

double  *geoTransform  ;//與前面的仿射變換的相似
GDALGCPsToGeoTransform( nGCPs, pGCPs, geoTransform, TRUE ); //由GCPs獲得仿射變換引數

int x,y;//任意點影像行列數

double xx,yy; //任意點影像對應的地理座標

double ex[nGCPs];//不同控制點算出的地理座標X
double ey[nGCPs];//不同控制點算出的地理座標Y

double pRadio[nGCPs]; //任意點相對與控制點的權值

double total=0; //權值和

double sx=0; //pRadio[i]*ex[i]
double sy=0; //pRadio[i]*ey[i]

for(int num=0;num<nGCPs;num++)
   {
        double gX,gY,tX,tY;
        tX=pGCPs[num].dfGCPPixel; //GCP位置X
        tY=pGCPs[num].dfGCPLine;  //GCP位置Y
        gX=pGCPs[num].dfGCPX;     //地理位置X
        gY=pGCPs[num].dfGCPY;     //地理位置Y

        pRadio[num]=1/(double)(sqrt((tX-x)*(tX-x)+(tY-y)*(tY-y)));
        total+=pRadio[num];

        ex[num]=gX+(x-tX)*geoTransform[1];
        ey[num]=gY+(y-tY)*geoTransform[5];

        sx+=pRadio[num]*ex[num];
        sy+=pRadio[num]*ey[num];
     }

xx=sx/total;

yy=sy/total;

對於遙感影象而言,如果想顯示在螢幕上,由於影象過大或者螢幕的關係,常常做成金字塔影像,或者拉伸影象資料,使影象產生變形以此來滿足顯示的需求。

假如已經將經緯度轉換成小數點位,即以度為單位,已知道螢幕的高(y)和寬(h),地理座標區域的範圍(maxLon,minLon,maxLat,minLat)..這裡我們知道了這些已知的引數

我們可以算出每畫素所代表的經度和緯度(有人稱這個為比例因子):
        公式:scaleX = h/((maxLon-minLon)*3600)  ----------X軸上每畫素代表的經度秒數;
        公式:scaleY = y/((maxLat-minLat)*3600)  -----------Y軸上每畫素代表的緯度秒數;

算出該地理座標區域中的任何一點(lon,lat)在螢幕上的座標

公式:screenX = lon*3600/scaleX;  ---------螢幕座標X軸座標
公式:screenY = lat*3600/scaleY; ----------螢幕座標Y軸座標,  

lon和lat為任意地點的地理座標,前面的仿射變換已經算出

假如需要佔滿整個螢幕:

 公式:minX = minLon*3600/scaleX;    區域左邊置最左端
公式:minY = minLat*3600/scaleY;     區域上面置最上端

當地地理範圍區域佔滿整個螢幕時,我們需要用到第三步計算出來的 screenX和screenY兩個引數
該區域中的任何一點的公式如下:   
 公式:X = screenX - minX = (lon - minLon)*3600/scaleX;     
由於緯度的方向和螢幕Y軸是相反的, 
 公式:screenMaxLat = (maxLat - minLat)*3600/scaleY;
 公式:screenLat = (lat - minLat)*3600/scaleY;
公式:Y = screenMaxLat - screenLat = (maxLat - lat)*3600/scaleY;