Proj
类主要是进行经纬度与地图投影坐标转换,以及反转。
可以参考前边对 proj
的介绍。 当初始化一个 Proj
类的实例时,地图投影的参数设置可以用关键字/值的形式。
关键字和值的形式也可以用字典或关键字参数,
或者一个 proj4 字符串(与 proj 命令兼容)。
如果可选的关键字“errcheck
”为真的话(默认为假),
一个异常将会被给出,如果转换无效的话。
如果为假,且转换无效时,没有异常抛出,会返会一个无效值 1.e30。
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))
proj4投影控制参数必须在字典projparams或关键字参数给定。看到项目文件(http://proj.maptools.org )有关指定投影参数的更多信息。
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)
lon,lat=p(x,y,inverse=True)
print('%.3f,%.3f' %(lon,lat))
import math #验证关键字 radians
x,y=p(math.radians(105),math.radians(36),radians=True)
print( '%.3f,%.3f' %(x,y) )
lons=(105,106,104)
lats=(36,35,34)
x,y=p(lons,lats) #将经纬度放入元组中
print('%.3f,%.3f,%.3f' %x) #普通打印
print('%.3f,%.3f,%.3f' %y)
type(x) #输出的类型为元组
zip(x,y) #用 zip函数包装
utm=Proj(proj='utm',zone=48,ellps='WGS84')
x,y=utm(105,36) #用关键字定义一个投影。
x,y
utm.is_geocent() #很纳闷,看来有必要弄清什么是 geocentric
utm.is_latlong()
latlong=Proj('+proj=latlong')
latlong.is_latlong()
latlong.is_geocent()
投影转换
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
都是Proj
对象。
在 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
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
utm_x,utm_y=utm(105,36)
print(utm_x,utm_y )
# 下边直接从 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 )
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)
实例:进行两个点的等分
例如:
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))
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.
az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
print ("%7.3f %6.3f %12.3f" % (az12,az21,dist))
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 %13.3f" % (endlat,endlon,backaz))
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))
import helper; helper.info()