Proj 类

Proj 类主要是进行经纬度与地图投影坐标转换,以及反转。 可以参考前边对 proj 的介绍。 当初始化一个 Proj 类的实例时,地图投影的参数设置可以用关键字/值的形式。 关键字和值的形式也可以用字典或关键字参数, 或者一个 proj4 字符串(与 proj 命令兼容)。

如果可选的关键字“errcheck”为真的话(默认为假), 一个异常将会被给出,如果转换无效的话。 如果为假,且转换无效时,没有异常抛出,会返会一个无效值 1.e30。

转换

首先看一下如何转换成平面坐标:

In [4]:
from pyproj import Proj
p=Proj('+proj=aea +lon_0=105 +lat_1=25 +lat_2=47 +ellps=krass')
x,y=p(105,36)
print('%.3f,%.3f' %(x,y))
0.000,3847866.973

proj4投影控制参数必须在字典projparams或关键字参数给定。看到项目文件(http://proj.maptools.org )有关指定投影参数的更多信息。

In [5]:
Proj(proj='utm',zone=10,ellps='WGS84') # use kwargs
Proj('+proj=utm +zone=10 +ellps=WGS84') # use proj4 string
Proj(init="epsg:32667")
Proj("+init=epsg:32667",preserve_units=True)
Out[5]:
<pyproj.Proj at 0x7fc7f27353a8>

当调用一个经纬度的 Proj 类的实例时, 将会把十进制的经纬度,转换成为地图的 (x, y) 投影。

上面第1行导入Proj类; 第2行则使用 proj4 格式初使化一个投影(中国等积投影); 第3行则是进行转换。

逆变换

In [6]:
lon,lat=p(x,y,inverse=True)
print('%.3f,%.3f' %(lon,lat))
105.000,36.000

如果可选的关键字 inverse 等于 True 的时候(默认为假),则进行相反的转换。

第1行重新进行转换; 第2行则是打印输出转换的结果。 从结果来看,与最开始的输入结果是一致的。

弧度变换

pyproj.Proj类除了使用度作为单位输入之外,pyproj.Proj类还接受弧度作为单位的输入。

In [7]:
import math   #验证关键字 radians
x,y=p(math.radians(105),math.radians(36),radians=True) 
print( '%.3f,%.3f' %(x,y) )
0.000,3847866.973

使用math模块,将度转换为弧度。

如果关键字 radiansTrue 的话(默认为假), 则经纬度的单位则是弧度,而不是度。

使用数组进行批量转换

pyproj.Proj类除进行单个值的处理, pyproj.Proj 类还接受数组的输入,进行批量转换。 这个在实际的开发中是非常实用的。

In [8]:
lons=(105,106,104)   
lats=(36,35,34) 
x,y=p(lons,lats)  #将经纬度放入元组中
print('%.3f,%.3f,%.3f' %x)    #普通打印
0.000,89660.498,-90797.784
In [9]:
print('%.3f,%.3f,%.3f' %y)
3847866.973,3735328.476,3622421.811
In [10]:
type(x)    #输出的类型为元组  
Out[10]:
tuple
In [11]:
zip(x,y)   #用 zip函数包装    
Out[11]:
<zip at 0x7fc7f2771808>

可以将经纬度分别存入一个 list或 array。可以进行更高效率的转换。 输入的值应当是双精度(如果输入的不是,它们将会被转为双精度)。

虽然 Proj 可以与numpy和常规python数组对象, python序列和标量,但是使用 array 对象速度快一些

使用关键字定义投影

pyproj除了使用 proj4 字符串进行投影的定义之外, pyproj还支持使用关键字来定义投影。

In [12]:
utm=Proj(proj='utm',zone=48,ellps='WGS84') 
x,y=utm(105,36)   #用关键字定义一个投影。
x,y 
Out[12]:
(499999.99999999773, 3983948.4533356656)

初始化一个投影: Proj4投影控制参数或者是以字典形式给出, 或者是以关键字参数给出,也可以用 proj4的形式给出字符串。

两个 Proj实例的函数

is_geocent(self) 返回 True当投影为 geocentric (x/y) coordinates。
is_latlong(self) 返回 True当为地理坐标系经纬度时。

举例:

In [13]:
utm.is_geocent() #很纳闷,看来有必要弄清什么是 geocentric  
Out[13]:
False
In [14]:
utm.is_latlong()  
Out[14]:
False
In [15]:
latlong=Proj('+proj=latlong') 
latlong.is_latlong() 
Out[15]:
True
In [16]:
latlong.is_geocent() 
Out[16]:
False

投影转换

transform()函数是在pyproj库下, 可以进行两个不同投影的转换。 相当于proj程序里的cs2cs子程序。用法如下:

transform(p1, p2, x, y, z=None, radians=False)
x2, y2, z2 = transform(p1, p2, x1, y1, z1, radians=False)

这里p1p2都是Proj对象。

p1p2两个投影之间进行投影转换。把 p1坐标系下的点(x1,y1,z1)转换到 p2所定义的投影中去。 z1是可选的,如果没有设 z1的话,将假定为 0, 并返回 x2y2。 使用这个函数的时候要注意不要进行基准面的变换(datum)。 关键字radians只有在 p1p2中有一个为地理坐标系时才起作用, 并且是把地理坐标系的投影的值当作弧度。 判断是否为地理坐标可以用 p1.is_latlong()p2.is_latlong()函数。输入的 (x,y,z)可以分别是数组或序列的某一种形式。

例如:

In [17]:
# import pyproj 
from pyproj import Proj
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 
Out[17]:
(0.0, 3847866.972516728)
In [18]:
utm_x,utm_y=utm(105,36) 
print(utm_x,utm_y )
499999.99999999773 3984019.058813517
In [19]:
# 下边直接从 albers转为 utm坐标
from pyproj import transform
to_utm_x,to_utm_y = transform(albers,utm,albers_x ,albers_y) 
print(to_utm_x,to_utm_y )
499999.99999999773 3984019.058813516

Geod类

Geod类主要用来求算地球大圆两点间的距离和相对方位,以及相反的操作。 同时也可以在两点间插入等分点。

该类主要包括三个函数:

正转换

fwd()函数可以进行正转换,返回经纬度。也可以用 NumPy 与常规的Python数组对象,Python序列和标量。如果 弧度为真,把输入的经纬度单位当作弧度,而不是度。距离单位为米。

fwd(self,lons,lats,az,dist,radians=False)

参数为经度,纬度,相对方位,距离。

逆变换

inv()函数可以进行逆变换,已知两点经纬度,返回其前方位解,后方位角,以及距离。

inv(self, lons1, lats1, lons2, lats2, radians=False)

弧线等分

npts()函数可以进行弧线等分。给出一个始点(lon1,lat1)和终点(lon2,lat2)以及等分点数目 npts。

npts(self, lon1, lat1, lon2, lat2, npts, radians=False)

实例:进行两个点的等分

例如:

In [20]:
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) #单位为米

这样对于一个确定的实例就可以使用上边的三个函数了。

例如:

In [21]:
from pyproj import Geod
g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.

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.

In [22]:
az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
print ("%7.3f %6.3f %12.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.

In [23]:
endlon, endlat, backaz = g.fwd(boston_lon,boston_lat, az12, dist)
print("%6.3f  %6.3f %13.3f" % (endlat,endlon,backaz))
45.517  -123.683        75.654

compute the azimuths, distances from New York to several cities (pass a list)

In [24]:
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.651

In [25]:
import helper; helper.info()
页面更新时间: 2019-02-28 18:41:06
操作系统/OS: Linux-4.9.0-8-amd64-x86_64-with-debian-9.8
Python: 3.5.3
In [ ]: