PROJ4初探(轉並整理格式)
Proj4是一個免費的GIS工具,軟件還稱不上。 它專注於地圖投影的表達,以及轉換。采用一種非常簡單明了的投影表達--PROJ4,比其它的投影定義簡單,但很明顯。很容易就能看到各種地理坐標系和地 圖投影的參數,同時它強大的投影轉換功能,也是非常吸引人的。許多的 GIS軟件中也將其集成在內。Proj可以在 window的命令下有可運行的 EXE文件,其實它更主要的是一個庫!可以用來編一些批處理。在 Linux下除了可以直接運行外,還可以作為庫來進行更高功能的開發。
1 安裝
Window下安裝:
從Proj4
的網站下載安裝文件,解壓縮,把路徑加到環境變量里即可。具體操作步驟在解壓后的README有詳細說明。
Ubuntu Linux下安裝: 可以在終端輸入:
sudo apt -get install proj
Fedora Linux 下安裝: 可以在終端輸入:
su -c 'yum install proj'
2 快速開始
在終端或 DOSshell下可以輸入(帶$的為向終端里輸入的命令):
$proj Rel. 4.6.0, 21 Dec 2007 usage: proj [-beEfiIlormsStTvVwW [args]] [+opts[=arg]] [files]
會顯示出 proj的用法。包括參數設置,可選項,和輸入文件。
2.1 顯示參數
對於作地圖和 GIS工作者來說投影可謂是一切的基礎,投影的正確與否將關第到最終結果正確與否。在 proj里邊集成了,許多的制作地圖用的投影參數。我們可以使用下邊的命令來顯示在 proj里的內置的有關地圖投影的參數。顯示投影類型:
$proj -l aea: Albers Equal Area aeqd: Azimuthal Equidistant ... ... ... wag5: Wagner V wag6: Wagner VI wag7: Wagner VII weren: Werenskiold I wink1: Winkel I wink2: Winkel II wintri: Winkel Tripel
同樣的,還有命令:
$proj -le
顯示支持的橢球體(ellipsoid)信息,顯示結果省略。
$proj -ld
顯示Proj4支持的基准面(Datum)信息,顯示結果省略。
這兩個概念是有區別的。學過地圖學的都知道,地圖學上對地球上的抽象,第一次抽象為水准面
(等重力面),第二次抽象為橢球體
(ellipsoid),第三次 抽象現在我認為是將橢球體進行定位之后,所確定的具有明確的方向的橢球體,它的要求能夠很好的為當地區的地圖制作服務,這個似乎才可稱為基准面
(Datum)。當輸入:
$ proj -ld __datum_id__ __ellipse___ __definition/comments______________________________ WGS84 WGS84 towgs84=0,0,0 GGRS87 GRS80 towgs84=-199.87,74.79,246.62 Greek_Geodetic_Reference_System_1987 NAD83 GRS80 towgs84=0,0,0 North_American_Datum_1983 NAD27 clrk66 nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat
上邊顯示的就是基准面和橢球體的差異。
2.2 投影轉換
我 國常用的地圖投影主要有:Albers,Lambert,Gauss-Kruger,UTM投影。 等積投影由於沒有面積變形,所以在土地調查,植被蓋度分類等涉及到要保持面積不能變形的情況.中國的全國性地圖許多采用等積投影。國際上稱為Albers
投影,是一種圓錐等積投影。中國所使用的 Albers的參數是雙標准緯線
,25N
,47N
,中央經線為105E
,橢球體為Krassovsky
。用proj4表示為:
+proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47
下邊將用中國的 Albers投影,簡稱為 Albers_China來作個簡單的投影轉換。
$proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 105 36 0.00 3847866.97 104d36'54 36d25'9 -33897.90 3895309.74 104d25'36.9E 36d52'41N -50158.40 3947261.73
也可以進行批量轉:
$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 <<EOF > 105 36 > 104 36 > 106 24 > EOF 0.00 3847866.97 -88522.43 3848312.80 102064.08 2503934.26
同樣也可以進行反轉,即將 Albers轉為經緯度,只要在命令中加入參數-I
$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 -I <<e > 0 3847866.97 > -88522.43 3848312.80 > 102064.08 2503934.26 > e 105dE 36dN 104dE 36dN 106dE 24dN
在這里轉換的過程中始終是按經度-緯度,x-y的順序放進的。你也許會想將它們的方向掉轉。如果是輸入時想轉可在命令中加 -r,如果是輸出想掉轉,可以是加 -s
$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 -r -s <<e > 36 105 > 33 104 > e 3847866.97 0.00 3509623.92 -91933.97
同樣也可以通過文件來進行批量轉換:
lat_lon.test 105dE 36dN 104dE 36dN 106dE 24dN $ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 ~/lat_lon.test >alberst.test
生成的:alberst.test
0.00 3847866.97 -88522.43 3848312.80 102064.08 2503934.26
你也可以在文件中加注釋和對坐標點的說明,在轉換后仍可以保留:
lat_lon.test #it's just a test for convert file format 105dE 36dN not Lanzhou 104dE 36dN Lanzhou 106dE 24dN Unknow place
命令:
$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 ~/lat_lon.test >albers.test albers.test
#it's just a test for convert file format 0.00 3847866.97 not Lanzhou -88522.43 3848312.80 Lanzhou 102064.08 2503934.26 Unknow place
在命令上邊的~/lat_lon.test
是輸入的文件,~
在 linux下指的是當前目錄,windows下沒試過,不過可以用絕對路徑。>
是重定向,輸出文件。
2.3 地圖單位
proj4支持許多單位,可以通過proj -lu
,看到支持的單位:
$ proj -lu km 1000. Kilometer m 1. Meter dm 1/10 Decimeter cm 1/100 Centimeter mm 1/1000 Millimeter kmi 1852.0 International Nautical Mile
其中 proj默認的單位為米
(meter),我們設置參數+units
來控制輸入的坐標單位。我們可以將它的輸入或輸出的數據的單位改為其它:
$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 +units=km -I <<e > 0 3847.86697 > -88.52243 3848.31280 > 102.06408 2503.93426 > e 105dE 36dN 104dE 36dN 106dE 24dN
3 地圖投影設置
其實在上邊我們已經用到了一些參數。比如進行投影的反轉,所使的是-I。還有就是給出的中國等積投影是:
proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25
其 中有許多參數前邊都加了前綴+
,這后邊的參數是對地圖投影的真正的設置,proj命令也將按這個規定來進行轉換。這種的參數的形式是+param=value
為 param(參數)來給定一個value(值),這個值可以是一個以度分秒格式或實數,整數,甚至可以是一個ASCII字符串。所要注意的是一種拼錯了的 參數名會被不理會,如果一個參數輸入了兩次,則第一次輸入的將被使用。所以在執行前要將參數檢查好。另一個很有用的特點是,proj會自動決定中央經線並 且追加一個+lon_0參數到你的定義中,如果你沒有定義中央經線的話。
3.1 選擇投影
選擇投影使用的參數是+proj=name
。proj中集成了不少的投影而且還在不斷的增加。可惜的對我們中國的並沒增加多少!它在這里提供我們增加我們自己投影的方法。我們可以用命令:
$proj -l #得到內置的投影類型
aea : Albers Equal Area #可以看到 Albers等積投影的在 proj中的值是 aea
$proj -lP #得到更詳細的投影信息,投影包含的參數。
aea : Albers Equal Area #更詳細的信息
Conic Sph&Ell lat_1= lat_2=
這樣我們就可以來定你想要的投影了。
3.2 選擇橢球體
雖然投影定了,可是還要確定橢球體。有兩種辦法:一種是利用內置的一些橢球體。
$proj -le #顯示出內置的橢球體。
找幾個我國經常用的橢球體:
clrk66 a=6378206.4 b=6356583.8 Clarke 1866 krass a=6378245.0 rf=298.3 Krassovsky, 1942 WGS84 a=6378137.0 rf=298.257223563 WGS 84
另一種就是你可以自定義一個你自己的地球,在上邊的三個橢球體上就包含了定義一個橢球體的參數。以下面是這幾個參數的含義:必須有的:
+a=A
橢球的赤道半徑(半長軸)
下邊的這些只要一個就行的參數:
+b=B
橢球的極半徑(半短軸)
+f=F
橢球的扁率F=(A-B)/A +rf=RF
橢球的反扁率RF=1/F +e=E
偏心率+es=ES
偏心率的平方E^2
例如指定一個 Clark1866 橢球來作為中國等積投影的參數就可以這樣來設置:
$proj +proj=aea +ellps=clrk66 +lon_0=105 +lat_1=25 105 36 #可以和上邊的 krass橢球體比較 0.00 3850517.66 104 36 -88731.89 3850957.52 106 24 102046.72 2507997.23
同樣也可以自定義:
$ proj +proj=aea +a=6378206.4 +es=.006768658 +lon_0=105 +lat_1=25
或:
$ proj +proj=aea +a=6378206.4 b=6356583.8 +lon_0=105 +lat_1=25
得到結果都是一樣的。如果你僅僅只指定了一個+a則會得到一個正規的球體。例如:
proj +a=1
你會定義一個半徑R=1
的單位球體。
3.3 通用的參數
我們知道UTM投影
,為確保每個帶中的點全為正,在北半球將坐標從每帶的中央徑線西移500
公里。而南半球為了保證Y軸為正,不得不向南移。這樣兩個參數在 proj中可以用兩個參數+x_0
和+y_0
來確定。例如在UTM中,
+x_0=5000000 +y_0=0
還有一個是+lat_0
:不太清楚,把原文貼過來。
A fourth parameter, lat_0= 0, is used to designate a φcentral parallel and associated y axis origin for several projections.
還有最后一個+lon_0
:指中央徑線,在 UTM還有 Albers中要設置。上邊的四個參數是可選的如果沒有將會默認為 0
,除+lon_0
,它會計算。這樣我們就完成了一個投影的完整定義了,至少我認為這是一種很很簡單很優美的定義投影的方式。通常把這種proj定義出的投影形式稱為PROJ4格式
。例如剛開始寫的中國的等積投影:
+proj=aea +ellps=krass +lon_0=105 +lat_1=25
還有中國的Lambert
投影,寫法和上邊差不多,在 proj4里的投影名稱叫+proj=lcc
定義一個 UTM 投影的方法
$proj +proj=utm +zone=48 105 36 #在離蘭州最近的一條中央經線 500000.00 3983948.45 #可以看出西移 500000米
定義一個高斯克呂格(它是橫軸墨卡托的一種,也是橫軸黑卡托的默認值)
$proj +proj=tmerc +lon_0=105 #這上邊沒帶,可以簡單的算出每個帶的中央經線的度數:(帶數*6)-3就是中央經線
4 自定義你的 proj 及內置的其它程序介紹
4.1 自定義你的參數
也許我們或許希望將我們常用的地圖投影的定義放到文件中,我們只要在需要時將其導入即可。Proj4也提供了這個功能。我們以設置+init
參數來進行。 首先先將地圖投影保存到文件中。我們保存地圖投影的文件格式如下: proj.dat文件
中國全國范圍內的Albers投影
<albers_china> +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 <>
上 邊我們將中國的等積投影保存在名為proj.dat
的文件中。作為這個放置投影的文件的格式是:#
表示注釋,proj
不執行 <...>
夾在這里邊的是關鍵字,因為在一個投影文件中可以放多個投影, proj在這個文件里靠關鍵字來選擇投影。<>
最的這個表示投影定義的結束。中間的部分就表示為定義的投影了。這樣我們就可以很簡單的使用+init
來使用我們定義的這個投影了!
$ proj +init=~/proj.dat:albers_china 106 36 88522.43 3848312.80 105 36 0.00 3847866.97 106 24 102064.08 2503934.26
+init
后邊是個等號,之后是你保存投影文件的絕對路徑。之后是個冒號,冒號之后則是投影的關鍵字。由於我是在 LINUX上的,~
號代表的是當前目錄。在 Windows下則要改為類似於:+init=c:\proj.dat:albers_china
其實在 proj初使時它已經有一些默認的值了,比如橢球體默認為WGS84
。如果要更改這些。可以參考后邊所列的參考書。
4.2 proj4 里的其它程序
其實在你安裝完 proj4之后,不光只安裝了個 proj程序,還包括其它的幾個,以及一個 lib。先介紹其內置的幾個其它函數。
invproj
反轉,將xy
轉換為lonla
t。相當於proj -I
cs2cs
是從一個投影轉到另一個投影,具體的使用方法是:
$ cs2cs +init=~/proj.dat:albers_china +to +proj=utm +zone=48 +ellps=krass 88522.43 3848312.80 590130.56 3984481.41 0.00 +to
前邊你輸入的坐標所屬的投影,后邊是要轉出的投影。這就解決了,如果要用proj
在兩個坐標系間轉,必須要先轉為經緯度。其它的參數可以參照proj的設置。
nad2nad
這是他們給他們自己國家的兩個基准面進行基准面轉換所作的。
轉個原文吧,不知道我們什么時候用上自己的基准面轉換,從54轉80
Program nad2nad is a filter to convert data between North America Datum 1927 (NAD27) and North American Datum 1983. nad2nad can option-ally process both State Plane Coordinate System (SPCS) and Universal Transverse Mercator (UTM) grid data as well as geographic data for both input and output. This can also be accomplished with the cs2cs program.
geod
這是個很有意思的,在定義了橢球體之后可以求出地表任意兩點的大圓距離,以及兩點的相對方位。當然也可以有它的反函數,invgeod
。其使用方法,可以看幫助,一目了然。
5 Proj 的參數全解
下邊將 Proj程序進行一個小的總結,對它所包含的所有參數進行說明。投影變換的命令為:
proj [-control] [+control] [files]
在 unix,Linux上還包括一個反變換:
invproj [files]
代表輸入的數據文件,執行順序是從左向右,如里這里為空的話,則程序默認是從標准輸入設備中輸入,即鍵盤。 [-control] 參數被用來制約控制,數據的輸入輸出,和選擇要計算的基本信息。下邊將每個-control
參數進行說明。
-I 進行投影的反轉換。在上邊說到的 unix下的有反轉程序 invproj在windows下只要加上它,就可以進行反轉。從投影坐標轉為地理坐標。
-l[p|P|e|u|d] 顯示出 proj內置的一些地球參數。
-lp顯示投影。-lP顯示投影里的參數比較詳細。
-le顯示內置的橢球體,看完之后,基本上夠用的。
-lu顯示 proj支持的長度單位。
-ld顯示基准面,很少基本上都是美國的,也只有 WGS84我們現在用的比較多。
-b 用來把從標准輸入輸出中的數據當做二進制。輸入數據被認為是雙精度浮點。當 proj作為一個子程序時非常有用,可以忽略格式操作。
-i 僅選擇輸入為二進制,參考(-b)
-o 僅選擇輸出為二進制,參考(-b)
-ta 指定一個字符來作為注釋來不被執行,僅只能使用 ASCII 字符。(默認為#號,例如可以改為-t%,這就將%作為注釋符)
-e string 來作為當輸出錯誤時輸出的信息。注意 -e 空一格,之后再輸入你想要達到的輸出時輸出的信息。默認的錯誤輸出是 *\t*。所要注意的是如果選了-b,-o,-i的話,錯誤的消息將被使用系統的 HUGE VAL 輸出。
-r 這個用來反轉輸入坐標的順序,如果使用的話會將輸入默認的longtitude-latitude或者x-y改變為latitude-lontitude 或y-x。
-s 將輸出的順序進行反轉。參考 -r
-m 是指將輸出的值按比例進行縮小放,主要用來用來進行單位轉換,因為proj 默認的輸出單位是米,如果你想輸出為千米時就可以用它來指定。例如:-m 1:1000 或使用 -m 1/1000 都是可以的。也可將其放大使用 -m 1:0.1 顯示為分米為最小單位。
-f format 來對輸出進行精度規范,輸出的小數點精度。默認的格式為%.2f,即輸出兩位小數。對於經緯度則采用十進制來表示。例如:將輸出的小點數位保留五位,-f%.5f。
-[w|W]n 用來控制當使用度分秒格式時,秒的位數。n 是個小於10的數。
-wn 指保留的秒的精度是n位,如果在某一位后全為0時就不輸出。
-Wn 則強制將后邊的0補上,要求小數位為一個定長。這個有時很有用,在規定格式方面。
-v 來顯示你的選擇的投影的信息
-V 來顯示更詳細的你選擇投影的信息,增加了橢球體。主要用來檢查你的輸入是否是你想定義的投影,保證沒有錯誤。
-E
-S 來把每一點在投影后的各種變形顯示出來。將這些信息顯示在<>里。分別包括 h,k,s,omega,a 和 b。
[+control] 有加號的地圖參數都和地圖的投影有關,並且不同的投影對其中的有些參數不包括。除了+init 外,其它的帶+號的參數都可以作為初始文件中,或默認的文件。這些選擇是從按從左到右的順序進行處理的。當出現重復時,以前邊的參數為准。
+proj=name 這是經常必選的一個地圖投影轉換函數。Name 是一個投影名稱。
+init=file:key 名為 file 的文件使用 key 關鍵字,選擇包含了控制參數的投影。
+R=R 當把地球當成球體計算時,所要設置的半徑。
+ellps=acronym 選擇一個 proj 里的橢球體。Acronym 為橢球體名。
+a=a 設置橢球體的半長軸為a,單位為米
+es=e^2 設置橢球體的偏心度的平方。同樣的也可以用+b=b,+e=e,
+rf=1/f,+f=f。來設置,橢球體參數。可以參考上邊的說明。
+x_0=x0 設置假東,例如在gauss-kruge 中,東偏 5000000 米。主要用來保持坐標的非負。
+y_0=y0 設置假北,主要用來保持坐標非負。
+lon_0=c 高置中央經線。通常和+lat_0,一起決定投影的地理起點。
+lat_0=d 參考上邊的。
+units=name 設置地圖坐標的單位。
+geoc 地理坐標將被認為是地心坐標,當使用這個參數的時候。
+over 使用該參數時,把經緯不限制在-180~+180 間。
舉例:
$ proj +proj=goode 389 36 2611721.44 4007501.67 29 36 #當不使用時,會自動將經度換到范圍 2611721.44 4007501.67
$ proj +proj=goode +over #使用關鍵字時,則不會轉換。 389 36 35033090.98 4007501.67 29 36 2611721.44 4007501.67
在安裝Proj4的位置,有一個默認的文件夾,來指示不滿足要求的初始化文件名稱,以及投影的默認文件proj_def.dat
。環境變量PROJ_LIB
可以用來給出另一個新的文件夾。
6 在Python里使用的Proj4
proj 不光是一些應用程序的集合,它更是一個庫,其它語言可以來調用它,來進行更高級的開發和應用。在 proj安裝上之后,它本身作為庫,可以被 C/C++來調用。而 proj本身是一個開源的項目,同時 Python也是一個開源的編程語言。Proj理所應當的能夠用在 Python里。在 Python里的 Proj庫稱為 Pyproj。在 windows和 Linux下都很好安裝。
Pyporj
是 Python下的 proj。可以很方便的對點來進行地圖投影轉換。同時在它的基礎上開發出更高級的應用。Pyproj包里包括兩個類,Proj
類和Geod
類。
Proj
類相當於前面所說的proj
的功能。可以進行地圖投影的變換從經緯度轉為 xy投影坐標,也可以反轉。也可以在不同的地圖投影之間轉換。
Geod
類相當於前邊介紹的proj
里的一個應用程序Geod
。可以很方便的計算地球上任意兩點的大圓距離,以及它們的相對方位。同時,也可根據方位和大圓距離來反算出另一點的經緯度。其處理的輸入坐標可以是 python數組,list\元組,scalar 或者 numpy/Numeric/numarray arrays。在導入 Pyproj后可以用其內部的函數 test()會運行一些例子。下邊介紹 Proj類和 Geod類,以及 transform函數。
6.1 Proj 類
Proj類主要是進行經緯度與地圖投影坐標轉換,以及反轉。可以參考前邊對 proj的介紹。當初始化一個 Proj類的時例時,地圖投影的參數設置可以用關鍵字\值的形式。關鍵字和值的形式也可以用字典,或關鍵字參數,或者一個 proj4 字符串(與 proj的命令兼容)。
當調用一個包含經緯度的 Proj類的實例時,將會把十進制的經緯度,轉換成為地圖的 xy坐標。如果可選的關鍵字inverse
等於 True的時候(默認為假),則進行相反的轉換。如果關鍵字'radians'為 True
的話(默認為假),則經緯度的單位則是弧度,而不是度。如果可選的關鍵字errcheck
為真的話(默認為假),一個異常將會被給出,如果轉換無效的話。如果為假的話,且轉換無效時,沒有異常拋出,但會返會一個無效值1.e30
。
可以將經緯度分別存入一個list或
array。可以進行更高效率的轉換。輸入的值應當是雙精度(如果輸入的不是,它們將會被轉為雙精度)。
雖然Proj
可以和numpy and regular python array objects,python sequences and scalars
,但是用array
對象速度快一些。
初始化一個投影:Proj4
投影控制參數或者是以字典形式給出,或者是以關鍵字參數給出,也可以用proj4
的形式給出字符串。例:
>>> from pyproj import Proj #從 pyproj導入 Proj類
>>> p=Proj('+proj=aea +lon_0=105 +lat_1=25 +lat_2=47 +ellps=krass') #初使化一個投影-中國等積投影,使用 proj4格式
>>> x,y=p(105,36) #進行轉換
>>> print '%.3f,%.3f' %(x,y) #按格式輸出 -0.000,3847866.973
>>> lon,lat=p(x,y,inverse=True) #進行投影反轉
>>> print '%.3,%.3' %(lon,lat) #xy為上邊轉的數值
>>> print '%.3f,%.3f' %(lon,lat) 105.000,36.000
>>> import math #驗證關鍵字 radians
>>> x,y=p(math.radians(105),math.radians(36),radians=True)
>>> print '%.3f,%.3f' %(x,y) -0.000,3847866.973
>>> lons=(105,106,104) >>> lats=(36,35,34)
>>> x,y=p(lons,lats) #將經緯度放入元組中
>>> print '%.3f,%.3f,%.3f' %x #普通打印 -0.000,89660.498,-90797.784
>>> print '%.3f,%.3f,%.3f' %y 3847866.973,3735328.476,3622421.811
>>> type(x) #輸出的類型為元組 <type 'tuple'>
>>> zip(x,y) 用 zip函數包裝 [(-1.1262207710106308e-09, 3847866.9725167281), (89660.498484069467, 3735328.4764740192), (-90797.783903947289, 3622421.8109651795)]
>>> utm=Proj(proj='utm',zone=48,ellps='WGS84')
>>> x,y=utm(105,36) #用關鍵字定義一個投影。
>>> x,y (500000.00000000116, 3983948.4533358263)
兩個 Proj實例的函數:
is_geocent(self)
返回True
當投影為geocentric(x/y) coordinates
。
is_latlong(self)
返回True
當為地理坐標系經緯度時。
舉例:
>>> utm.is_geocent() #很納悶,看來有必要弄清什么是 geocentric False
>>> utm.is_latlong() False
>>> latlong=Proj('+proj=latlong')
>>> latlong.is_latlong() True
>>> latlong.is_geocent() False
6.2 transform()函數
transform()
函數是在pyproj
庫下,可以進行兩個不同投影的轉換。相當於proj
程序里的cs2cs
子程序。用法如下:
transform(p1, p2, x, y, z=None, radians=False) x2, y2, z2 = transform(p1, p2, x1, y1, z1, radians=False)
在 p1,p2兩個投影間進行投影轉換。將把在 p1坐標系下的點(x1,y1,z1)
轉換到 p2所定義的投影中去。z1
是可選的,如果沒有設 z1的話將會假定為 0,並僅僅返回 x2,y2
。使用這個函數的時候要注意不要進行基准面的變換(datum)。關鍵字radians
只有在 p1,p2中有一個為地理坐標系時才起作用,並且把是地理坐標系的投影的值會當作弧度。判斷是否為地理坐標可以用 p1.is_latlong()
和p2.is_latlong()
函數。 輸入的x,y,z
可以是分別是數組或序列的某一種形式。舉例:
>>> import pyproj
>>> albers=Proj('+proj=aea +lon_0=105 +lat_1=25 +lat_2=47 +ellps=krass')
>>> utm=Proj(proj='utm',zone=48,ellps='krass')
>>> albers_x,albers_y=albers(105,36)
>>> albers_x,albers_y (-1.1262207710106308e-09, 3847866.9725167281)
>>> utm_x,utm_y=utm(105,36)
>>> utm_x,utm_y (500000.00000000116, 3984019.0588136762)
#下邊直接從 albers轉為 utm坐標
>>>to_utm_x,to_utm_y=pyproj.transform(albers,utm,albers_x ,albers_y)
>>> to_utm_x,to_utm_y (500000.00000000116, 3984019.0588136753)
6.3 Geod類
主要用來求算地球大圓兩點間的距離,及其相對方位,以及相反的操作。同時也可以在兩點間插入 N等分點。
該類主要包括三個函數:
fwd(self,lons,lats,az,dist,radians=False)
參數為經度,緯度,相對方位,距離。
正轉換,返回經緯度,以及 back azimuths。可以用 numpy and regular python array objects, python
sequences and scalars。如果 radians為真把輸入的經緯度單位當作弧度,而不是度。距離單位為米。
inv(self, lons1, lats1, lons2, lats2, radians=False)
反變換,已知兩點經緯度,返回其前方位解,后方位角,以及距離。
npts(self, lon1, lat1, lon2, lat2, npts, radians=False)
給出一個始點(lon1,lat1)和終點(lon2,lat2)以及等份點數目 npts。舉例(懶得寫了,照抄過來,已驗證):
>>> from pyproj import Geod
>>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
>>> # specify the lat/lons of Boston and Portland.
>>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
>>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
>>> # find ten equally spaced points between Boston and Portland.
>>> lonlats = g.npts(boston_lon,boston_lat,portland_lon,portland_lat,10)
>>> for lon,lat in lonlats: print '%6.3f %7.3f' % (lat, lon) 43.528 -75.414 44.637 -79.883 45.565 -84.512 46.299 -89.279 46.830 -94.156 47.149 -99.112 47.251 -104.106 47.136 -109.100 46.805 -114.051 46.262 -118.924
初始化一個 Geod類實例:使用關鍵字方法來傳一個橢球體,橢球體為 proj中支持的任何一個。
>>> g=Geod(ellps='krass')
同時也可以指定 a,b,f,rf,e,es來設定地球橢球體。
>>> miniearth=Geod(a=2,b=1.97) #單位為米
這樣對於一個確定的實例就可以使用上邊的三個函數了。實例(繼續照抄過來):
>>> from pyproj import Geod
>>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
>>> # specify the lat/lons of some cities.
>>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
>>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
>>> newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)
>>> london_lat = 51.+(32./60.); london_lon = -(5./60.)
>>> # compute forward and back azimuths, plus distance
>>> # between Boston and Portland. >>> az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
>>> print "%7.3f %6.3f .3f" % (az12,az21,dist) -66.531 75.654 4164192.708 >>> # compute latitude, longitude and back azimuth of Portland, >>> # given Boston lat/lon, forward azimuth and distance to Portland. >>> endlon, endlat, backaz = g.fwd(boston_lon,boston_lat, az12, dist) >>> print "%6.3f %6.3f .3f" % (endlat,endlon,backaz) 45.517 -123.683 75.654 >>> # compute the azimuths, distances from New York to several >>> # cities (pass a list) >>> lons1 = 3*[newyork_lon]; lats1 = 3*[newyork_lat] >>> lons2 = [boston_lon, portland_lon, london_lon] >>> lats2 = [boston_lat, portland_lat, london_lat] >>> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2) >>> for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f %9.3f" % (faz,baz,d) 54.663 -123.448 288303.720 -65.463 79.342 4013037.318 51.254 -71.576 5579916.649
7 參考資料